Coverage for /builds/kinetik161/ase/ase/calculators/castep.py: 41.82%
1442 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +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
31from typing import List, Set
33import numpy as np
35import ase
36import ase.units as units
37from ase.calculators.calculator import (PropertyNotImplementedError,
38 compare_atoms, kpts2sizeandoffsets)
39from ase.calculators.general import Calculator
40from ase.config import cfg
41from ase.constraints import FixCartesian
42from ase.dft.kpoints import BandPath
43from ase.io.castep import read_bands, read_param
44from ase.parallel import paropen
46__all__ = [
47 'Castep',
48 'CastepCell',
49 'CastepParam',
50 'create_castep_keywords']
52contact_email = 'simon.rittmeyer@tum.de'
54# A convenient table to avoid the previously used "eval"
55_tf_table = {
56 '': True, # Just the keyword is equivalent to True
57 'True': True,
58 'False': False}
61def _self_getter(getf):
62 # A decorator that makes it so that if no 'atoms' argument is passed to a
63 # getter function, self.atoms is used instead
65 def decor_getf(self, atoms=None, *args, **kwargs):
67 if atoms is None:
68 atoms = self.atoms
70 return getf(self, atoms, *args, **kwargs)
72 return decor_getf
75def _parse_tss_block(value, scaled=False):
76 # Parse the assigned value for a Transition State Search structure block
77 is_atoms = isinstance(value, ase.atoms.Atoms)
78 try:
79 is_strlist = all(map(lambda x: isinstance(x, str), value))
80 except TypeError:
81 is_strlist = False
83 if not is_atoms:
84 if not is_strlist:
85 # Invalid!
86 raise TypeError('castep.cell.positions_abs/frac_intermediate/'
87 'product expects Atoms object or list of strings')
89 # First line must be Angstroms, or nothing
90 has_units = len(value[0].strip().split()) == 1
91 if (not scaled) and has_units and value[0].strip() != 'ang':
92 raise RuntimeError('Only ang units currently supported in castep.'
93 'cell.positions_abs_intermediate/product')
94 return '\n'.join(map(str.strip, value))
95 else:
96 text_block = '' if scaled else 'ang\n'
97 positions = (value.get_scaled_positions() if scaled else
98 value.get_positions())
99 symbols = value.get_chemical_symbols()
100 for s, p in zip(symbols, positions):
101 text_block += ' {} {:.3f} {:.3f} {:.3f}\n'.format(s, *p)
103 return text_block
106class Castep(Calculator):
107 r"""
108CASTEP Interface Documentation
111Introduction
112============
114CASTEP_ [1]_ W_ is a software package which uses density functional theory to
115provide a good atomic-level description of all manner of materials and
116molecules. CASTEP can give information about total energies, forces and
117stresses on an atomic system, as well as calculating optimum geometries, band
118structures, optical spectra, phonon spectra and much more. It can also perform
119molecular dynamics simulations.
121The CASTEP calculator interface class offers intuitive access to all CASTEP
122settings and most results. All CASTEP specific settings are accessible via
123attribute access (*i.e*. ``calc.param.keyword = ...`` or
124``calc.cell.keyword = ...``)
127Getting Started:
128================
130Set the environment variables appropriately for your system::
132 export CASTEP_COMMAND=' ... '
133 export CASTEP_PP_PATH=' ... '
135Note: alternatively to CASTEP_PP_PATH one can set PSPOT_DIR
136as CASTEP consults this by default, i.e.::
138 export PSPOT_DIR=' ... '
141Running the Calculator
142======================
144The default initialization command for the CASTEP calculator is
146.. class:: Castep(directory='CASTEP', label='castep')
148To do a minimal run one only needs to set atoms, this will use all
149default settings of CASTEP, meaning LDA, singlepoint, etc..
151With a generated *castep_keywords.json* in place all options are accessible
152by inspection, *i.e.* tab-completion. This works best when using ``ipython``.
153All options can be accessed via ``calc.param.<TAB>`` or ``calc.cell.<TAB>``
154and documentation is printed with ``calc.param.<keyword> ?`` or
155``calc.cell.<keyword> ?``. All options can also be set directly
156using ``calc.keyword = ...`` or ``calc.KEYWORD = ...`` or even
157``ialc.KeYwOrD`` or directly as named arguments in the call to the constructor
158(*e.g.* ``Castep(task='GeometryOptimization')``).
159If using this calculator on a machine without CASTEP, one might choose to copy
160a *castep_keywords.json* file generated elsewhere in order to access this
161feature: the file will be used if located in the working directory,
162*$HOME/.ase/* or *ase/ase/calculators/* within the ASE library. The file should
163be generated the first time it is needed, but you can generate a new keywords
164file in the currect directory with ``python -m ase.calculators.castep``.
166All options that go into the ``.param`` file are held in an ``CastepParam``
167instance, while all options that go into the ``.cell`` file and don't belong
168to the atoms object are held in an ``CastepCell`` instance. Each instance can
169be created individually and can be added to calculators by attribute
170assignment, *i.e.* ``calc.param = param`` or ``calc.cell = cell``.
172All internal variables of the calculator start with an underscore (_).
173All cell attributes that clearly belong into the atoms object are blocked.
174Setting ``calc.atoms_attribute`` (*e.g.* ``= positions``) is sent directly to
175the atoms object.
178Arguments:
179==========
181========================= ====================================================
182Keyword Description
183========================= ====================================================
184``directory`` The relative path where all input and output files
185 will be placed. If this does not exist, it will be
186 created. Existing directories will be moved to
187 directory-TIMESTAMP unless self._rename_existing_dir
188 is set to false.
190``label`` The prefix of .param, .cell, .castep, etc. files.
192``castep_command`` Command to run castep. Can also be set via the bash
193 environment variable ``CASTEP_COMMAND``. If none is
194 given or found, will default to ``castep``
196``check_castep_version`` Boolean whether to check if the installed castep
197 version matches the version from which the available
198 options were deduced. Defaults to ``False``.
200``castep_pp_path`` The path where the pseudopotentials are stored. Can
201 also be set via the bash environment variables
202 ``PSPOT_DIR`` (preferred) and ``CASTEP_PP_PATH``.
203 Will default to the current working directory if
204 none is given or found. Note that pseudopotentials
205 may be generated on-the-fly if they are not found.
207``find_pspots`` Boolean whether to search for pseudopotentials in
208 ``<castep_pp_path>`` or not. If activated, files in
209 this directory will be checked for typical names. If
210 files are not found, they will be generated on the
211 fly, depending on the ``_build_missing_pspots``
212 value. A RuntimeError will be raised in case
213 multiple files per element are found. Defaults to
214 ``False``.
215``keyword_tolerance`` Integer to indicate the level of tolerance to apply
216 validation of any parameters set in the CastepCell
217 or CastepParam objects against the ones found in
218 castep_keywords. Levels are as following:
220 0 = no tolerance, keywords not found in
221 castep_keywords will raise an exception
223 1 = keywords not found will be accepted but produce
224 a warning (default)
226 2 = keywords not found will be accepted silently
228 3 = no attempt is made to look for
229 castep_keywords.json at all
230``castep_keywords`` Can be used to pass a CastepKeywords object that is
231 then used with no attempt to actually load a
232 castep_keywords.json file. Most useful for debugging
233 and testing purposes.
235========================= ====================================================
238Additional Settings
239===================
241========================= ====================================================
242Internal Setting Description
243========================= ====================================================
244``_castep_command`` (``=castep``): the actual shell command used to
245 call CASTEP.
247``_check_checkfile`` (``=True``): this makes write_param() only
248 write a continue or reuse statement if the
249 addressed .check or .castep_bin file exists in the
250 directory.
252``_copy_pspots`` (``=False``): if set to True the calculator will
253 actually copy the needed pseudo-potential (\*.usp)
254 file, usually it will only create symlinks.
256``_link_pspots`` (``=True``): if set to True the calculator will
257 actually will create symlinks to the needed pseudo
258 potentials. Set this option (and ``_copy_pspots``)
259 to False if you rather want to access your pseudo
260 potentials using the PSPOT_DIR environment variable
261 that is read by CASTEP.
262 *Note:* This option has no effect if ``copy_pspots``
263 is True..
265``_build_missing_pspots`` (``=True``): if set to True, castep will generate
266 missing pseudopotentials on the fly. If not, a
267 RuntimeError will be raised if not all files were
268 found.
270``_export_settings`` (``=True``): if this is set to
271 True, all calculator internal settings shown here
272 will be included in the .param in a comment line (#)
273 and can be read again by merge_param. merge_param
274 can be forced to ignore this directive using the
275 optional argument ``ignore_internal_keys=True``.
277``_force_write`` (``=True``): this controls wether the \*cell and
278 \*param will be overwritten.
280``_prepare_input_only`` (``=False``): If set to True, the calculator will
281 create \*cell und \*param file but not
282 start the calculation itself.
283 If this is used to prepare jobs locally
284 and run on a remote cluster it is recommended
285 to set ``_copy_pspots = True``.
287``_castep_pp_path`` (``='.'``) : the place where the calculator
288 will look for pseudo-potential files.
290``_find_pspots`` (``=False``): if set to True, the calculator will
291 try to find the respective pseudopotentials from
292 <_castep_pp_path>. As long as there are no multiple
293 files per element in this directory, the auto-detect
294 feature should be very robust. Raises a RuntimeError
295 if required files are not unique (multiple files per
296 element). Non existing pseudopotentials will be
297 generated, though this could be dangerous.
299``_rename_existing_dir`` (``=True``) : when using a new instance
300 of the calculator, this will move directories out of
301 the way that would be overwritten otherwise,
302 appending a date string.
304``_set_atoms`` (``=False``) : setting this to True will overwrite
305 any atoms object previously attached to the
306 calculator when reading a \.castep file. By de-
307 fault, the read() function will only create a new
308 atoms object if none has been attached and other-
309 wise try to assign forces etc. based on the atom's
310 positions. ``_set_atoms=True`` could be necessary
311 if one uses CASTEP's internal geometry optimization
312 (``calc.param.task='GeometryOptimization'``)
313 because then the positions get out of sync.
314 *Warning*: this option is generally not recommended
315 unless one knows one really needs it. There should
316 never be any need, if CASTEP is used as a
317 single-point calculator.
319``_track_output`` (``=False``) : if set to true, the interface
320 will append a number to the label on all input
321 and output files, where n is the number of calls
322 to this instance. *Warning*: this setting may con-
323 sume a lot more disk space because of the additio-
324 nal \*check files.
326``_try_reuse`` (``=_track_output``) : when setting this, the in-
327 terface will try to fetch the reuse file from the
328 previous run even if _track_output is True. By de-
329 fault it is equal to _track_output, but may be
330 overridden.
332 Since this behavior may not always be desirable for
333 single-point calculations. Regular reuse for *e.g.*
334 a geometry-optimization can be achieved by setting
335 ``calc.param.reuse = True``.
337``_pedantic`` (``=False``) if set to true, the calculator will
338 inform about settings probably wasting a lot of CPU
339 time or causing numerical inconsistencies.
341========================= ====================================================
343Special features:
344=================
347``.dryrun_ok()``
348 Runs ``castep_command seed -dryrun`` in a temporary directory return True if
349 all variables initialized ok. This is a fast way to catch errors in the
350 input. Afterwards _kpoints_used is set.
352``.merge_param()``
353 Takes a filename or filehandler of a .param file or CastepParam instance and
354 merges it into the current calculator instance, overwriting current settings
356``.keyword.clear()``
357 Can be used on any option like ``calc.param.keyword.clear()`` or
358 ``calc.cell.keyword.clear()`` to return to the CASTEP default.
360``.initialize()``
361 Creates all needed input in the ``_directory``. This can then copied to and
362 run in a place without ASE or even python.
364``.set_pspot('<library>')``
365 This automatically sets the pseudo-potential for all present species to
366 ``<Species>_<library>.usp``. Make sure that ``_castep_pp_path`` is set
367 correctly. Note that there is no check, if the file actually exists. If it
368 doesn't castep will crash! You may want to use ``find_pspots()`` instead.
370``.find_pspots(pspot=<library>, suffix=<suffix>)``
371 This automatically searches for pseudopotentials of type
372 ``<Species>_<library>.<suffix>`` or ``<Species>-<library>.<suffix>`` in
373 ``castep_pp_path` (make sure this is set correctly). Note that ``<Species>``
374 will be searched for case insensitive. Regular expressions are accepted, and
375 arguments ``'*'`` will be regarded as bash-like wildcards. Defaults are any
376 ``<library>`` and any ``<suffix>`` from ``['usp', 'UPF', 'recpot']``. If you
377 have well-organized folders with pseudopotentials of one kind, this should
378 work with the defaults.
380``print(calc)``
381 Prints a short summary of the calculator settings and atoms.
383``ase.io.castep.read_seed('path-to/seed')``
384 Given you have a combination of seed.{param,cell,castep} this will return an
385 atoms object with the last ionic positions in the .castep file and all other
386 settings parsed from the .cell and .param file. If no .castep file is found
387 the positions are taken from the .cell file. The output directory will be
388 set to the same directory, only the label is preceded by 'copy_of\_' to
389 avoid overwriting.
391``.set_kpts(kpoints)``
392 This is equivalent to initialising the calculator with
393 ``calc = Castep(kpts=kpoints)``. ``kpoints`` can be specified in many
394 convenient forms: simple Monkhorst-Pack grids can be specified e.g.
395 ``(2, 2, 3)`` or ``'2 2 3'``; lists of specific weighted k-points can be
396 given in reciprocal lattice coordinates e.g.
397 ``[[0, 0, 0, 0.25], [0.25, 0.25, 0.25, 0.75]]``; a dictionary syntax is
398 available for more complex requirements e.g.
399 ``{'size': (2, 2, 2), 'gamma': True}`` will give a Gamma-centered 2x2x2 M-P
400 grid, ``{'density': 10, 'gamma': False, 'even': False}`` will give a mesh
401 with density of at least 10 Ang (based on the unit cell of currently-attached
402 atoms) with an odd number of points in each direction and avoiding the Gamma
403 point.
405``.set_bandpath(bandpath)``
406 This is equivalent to initialialising the calculator with
407 ``calc=Castep(bandpath=bandpath)`` and may be set simultaneously with *kpts*.
408 It allows an electronic band structure path to be set up using ASE BandPath
409 objects. This enables a band structure calculation to be set up conveniently
410 using e.g. calc.set_bandpath(atoms.cell.bandpath().interpolate(npoints=200))
412``.band_structure(bandfile=None)``
413 Read a band structure from _seedname.bands_ file. This returns an ase
414 BandStructure object which may be plotted with e.g.
415 ``calc.band_structure().plot()``
417Notes/Issues:
418==============
420* Currently *only* the FixAtoms *constraint* is fully supported for reading and
421 writing. There is some experimental support for the FixCartesian constraint.
423* There is no support for the CASTEP *unit system*. Units of eV and Angstrom
424 are used throughout. In particular when converting total energies from
425 different calculators, one should check that the same CODATA_ version is
426 used for constants and conversion factors, respectively.
428.. _CASTEP: http://www.castep.org/
430.. _W: https://en.wikipedia.org/wiki/CASTEP
432.. _CODATA: https://physics.nist.gov/cuu/Constants/index.html
434.. [1] S. J. Clark, M. D. Segall, C. J. Pickard, P. J. Hasnip, M. J. Probert,
435 K. Refson, M. C. Payne Zeitschrift für Kristallographie 220(5-6)
436 pp.567- 570 (2005) PDF_.
438.. _PDF: http://www.tcm.phy.cam.ac.uk/castep/papers/ZKristallogr_2005.pdf
441End CASTEP Interface Documentation
442 """
444 # Class attributes !
445 # keys set through atoms object
446 atoms_keys = [
447 'charges',
448 'ionic_constraints',
449 'lattice_abs',
450 'lattice_cart',
451 'positions_abs',
452 'positions_abs_final',
453 'positions_abs_intermediate',
454 'positions_frac',
455 'positions_frac_final',
456 'positions_frac_intermediate']
458 atoms_obj_keys = [
459 'dipole',
460 'energy_free',
461 'energy_zero',
462 'fermi',
463 'forces',
464 'nbands',
465 'positions',
466 'stress',
467 'pressure']
469 internal_keys = [
470 '_castep_command',
471 '_check_checkfile',
472 '_copy_pspots',
473 '_link_pspots',
474 '_find_pspots',
475 '_build_missing_pspots',
476 '_directory',
477 '_export_settings',
478 '_force_write',
479 '_label',
480 '_prepare_input_only',
481 '_castep_pp_path',
482 '_rename_existing_dir',
483 '_set_atoms',
484 '_track_output',
485 '_try_reuse',
486 '_pedantic']
488 def __init__(self, directory='CASTEP', label='castep',
489 castep_command=None, check_castep_version=False,
490 castep_pp_path=None, find_pspots=False, keyword_tolerance=1,
491 castep_keywords=None, **kwargs):
493 self.__name__ = 'Castep'
495 # initialize the ase.calculators.general calculator
496 Calculator.__init__(self)
498 from ase.io.castep import write_cell
499 self._write_cell = write_cell
501 if castep_keywords is None:
502 castep_keywords = CastepKeywords(make_param_dict(),
503 make_cell_dict(),
504 [],
505 [],
506 0)
507 if keyword_tolerance < 3:
508 try:
509 castep_keywords = import_castep_keywords(castep_command)
510 except CastepVersionError as e:
511 if keyword_tolerance == 0:
512 raise e
513 else:
514 warnings.warn(str(e))
516 self._kw_tol = keyword_tolerance
517 keyword_tolerance = max(keyword_tolerance, 2) # 3 not accepted below
518 self.param = CastepParam(castep_keywords,
519 keyword_tolerance=keyword_tolerance)
520 self.cell = CastepCell(castep_keywords,
521 keyword_tolerance=keyword_tolerance)
523 ###################################
524 # Calculator state variables #
525 ###################################
526 self._calls = 0
527 self._castep_version = castep_keywords.castep_version
529 # collects warning from .castep files
530 self._warnings = []
531 # collects content from *.err file
532 self._error = None
533 # warnings raised by the ASE interface
534 self._interface_warnings = []
536 # store to check if recalculation is necessary
537 self._old_atoms = None
538 self._old_cell = None
539 self._old_param = None
541 ###################################
542 # Internal keys #
543 # Allow to tweak the behavior #
544 ###################################
545 self._opt = {}
546 self._castep_command = get_castep_command(castep_command)
547 self._castep_pp_path = get_castep_pp_path(castep_pp_path)
548 self._check_checkfile = True
549 self._copy_pspots = False
550 self._link_pspots = True
551 self._find_pspots = find_pspots
552 self._build_missing_pspots = True
553 self._directory = os.path.abspath(directory)
554 self._export_settings = True
555 self._force_write = True
556 self._label = label
557 self._prepare_input_only = False
558 self._rename_existing_dir = True
559 self._set_atoms = False
560 self._track_output = False
561 self._try_reuse = False
563 # turn off the pedantic user warnings
564 self._pedantic = False
566 # will be set on during runtime
567 self._seed = None
569 ###################################
570 # (Physical) result variables #
571 ###################################
572 self.atoms = None
573 # initialize result variables
574 self._forces = None
575 self._energy_total = None
576 self._energy_free = None
577 self._energy_0K = None
578 self._energy_total_corr = None
579 self._eigenvalues = None
580 self._efermi = None
581 self._ibz_kpts = None
582 self._ibz_weights = None
583 self._band_structure = None
585 # dispersion corrections
586 self._dispcorr_energy_total = None
587 self._dispcorr_energy_free = None
588 self._dispcorr_energy_0K = None
590 # spins and hirshfeld volumes
591 self._spins = None
592 self._hirsh_volrat = None
594 # Mulliken charges
595 self._mulliken_charges = None
596 # Hirshfeld charges
597 self._hirshfeld_charges = None
599 self._number_of_cell_constraints = None
600 self._output_verbosity = None
601 self._stress = None
602 self._pressure = None
603 self._unit_cell = None
604 self._kpoints = None
606 # pointers to other files used at runtime
607 self._check_file = None
608 self._castep_bin_file = None
610 # plane wave cutoff energy (may be derived during PP generation)
611 self._cut_off_energy = None
613 # runtime information
614 self._total_time = None
615 self._peak_memory = None
617 # check version of CASTEP options module against current one
618 if check_castep_version:
619 local_castep_version = get_castep_version(self._castep_command)
620 if not hasattr(self, '_castep_version'):
621 warnings.warn('No castep version found')
622 return
623 if not local_castep_version == self._castep_version:
624 warnings.warn(
625 'The options module was generated from version %s '
626 'while your are currently using CASTEP version %s' %
627 (self._castep_version,
628 get_castep_version(self._castep_command)))
629 self._castep_version = local_castep_version
631 # processes optional arguments in kw style
632 for keyword, value in kwargs.items():
633 # first fetch special keywords issued by ASE CLI
634 if keyword == 'kpts':
635 self.set_kpts(value)
636 elif keyword == 'bandpath':
637 self.set_bandpath(value)
638 elif keyword == 'xc':
639 self.xc_functional = value
640 elif keyword == 'ecut':
641 self.cut_off_energy = value
642 else: # the general case
643 self.__setattr__(keyword, value)
645 def band_structure(self, bandfile=None):
646 from ase.spectrum.band_structure import BandStructure
648 if bandfile is None:
649 bandfile = os.path.join(self._directory, self._seed) + '.bands'
651 if not os.path.exists(bandfile):
652 raise ValueError(f'Cannot find band file "{bandfile}".')
654 kpts, weights, eigenvalues, efermi = read_bands(bandfile)
656 # Get definitions of high-symmetry points
657 special_points = self.atoms.cell.bandpath(npoints=0).special_points
658 bandpath = BandPath(self.atoms.cell,
659 kpts=kpts,
660 special_points=special_points)
661 return BandStructure(bandpath, eigenvalues, reference=efermi)
663 def set_bandpath(self, bandpath):
664 """Set a band structure path from ase.dft.kpoints.BandPath object
666 This will set the bs_kpoint_list block with a set of specific points
667 determined in ASE. bs_kpoint_spacing will not be used; to modify the
668 number of points, consider using e.g. bandpath.resample(density=20) to
669 obtain a new dense path.
671 Args:
672 bandpath (:obj:`ase.dft.kpoints.BandPath` or None):
673 Set to None to remove list of band structure points. Otherwise,
674 sampling will follow BandPath parameters.
676 """
678 def clear_bs_keywords():
679 bs_keywords = product({'bs_kpoint', 'bs_kpoints'},
680 {'path', 'path_spacing',
681 'list',
682 'mp_grid', 'mp_spacing', 'mp_offset'})
683 for bs_tag in bs_keywords:
684 setattr(self.cell, '_'.join(bs_tag), None)
686 if bandpath is None:
687 clear_bs_keywords()
688 elif isinstance(bandpath, BandPath):
689 clear_bs_keywords()
690 self.cell.bs_kpoint_list = [' '.join(map(str, row))
691 for row in bandpath.kpts]
692 else:
693 raise TypeError('Band structure path must be an '
694 'ase.dft.kpoint.BandPath object')
696 def set_kpts(self, kpts):
697 """Set k-point mesh/path using a str, tuple or ASE features
699 Args:
700 kpts (None, tuple, str, dict):
702 This method will set the CASTEP parameters kpoints_mp_grid,
703 kpoints_mp_offset and kpoints_mp_spacing as appropriate. Unused
704 parameters will be set to None (i.e. not included in input files.)
706 If kpts=None, all these parameters are set as unused.
708 The simplest useful case is to give a 3-tuple of integers specifying
709 a Monkhorst-Pack grid. This may also be formatted as a string separated
710 by spaces; this is the format used internally before writing to the
711 input files.
713 A more powerful set of features is available when using a python
714 dictionary with the following allowed keys:
716 - 'size' (3-tuple of int) mesh of mesh dimensions
717 - 'density' (float) for BZ sampling density in points per recip. Ang
718 ( kpoint_mp_spacing = 1 / (2pi * density) ). An explicit MP mesh will
719 be set to allow for rounding/centering.
720 - 'spacing' (float) for BZ sampling density for maximum space between
721 sample points in reciprocal space. This is numerically equivalent to
722 the inbuilt ``calc.cell.kpoint_mp_spacing``, but will be converted to
723 'density' to allow for rounding/centering.
724 - 'even' (bool) to round each direction up to the nearest even number;
725 set False for odd numbers, leave as None for no odd/even rounding.
726 - 'gamma' (bool) to offset the Monkhorst-Pack grid to include
727 (0, 0, 0); set False to offset each direction avoiding 0.
728 """
730 def clear_mp_keywords():
731 mp_keywords = product({'kpoint', 'kpoints'},
732 {'mp_grid', 'mp_offset',
733 'mp_spacing', 'list'})
734 for kp_tag in mp_keywords:
735 setattr(self.cell, '_'.join(kp_tag), None)
737 # Case 1: Clear parameters with set_kpts(None)
738 if kpts is None:
739 clear_mp_keywords()
741 # Case 2: list of explicit k-points with weights
742 # e.g. [[ 0, 0, 0, 0.125],
743 # [ 0, -0.5, 0, 0.375],
744 # [-0.5, 0, -0.5, 0.375],
745 # [-0.5, -0.5, -0.5, 0.125]]
747 elif (isinstance(kpts, (tuple, list))
748 and isinstance(kpts[0], (tuple, list))):
750 if not all(map((lambda row: len(row) == 4), kpts)):
751 raise ValueError(
752 'In explicit kpt list each row should have 4 elements')
754 clear_mp_keywords()
755 self.cell.kpoint_list = [' '.join(map(str, row)) for row in kpts]
757 # Case 3: list of explicit kpts formatted as list of str
758 # i.e. the internal format of calc.kpoint_list split on \n
759 # e.g. ['0 0 0 0.125', '0 -0.5 0 0.375', '-0.5 0 -0.5 0.375']
760 elif isinstance(kpts, (tuple, list)) and isinstance(kpts[0], str):
762 if not all(map((lambda row: len(row.split()) == 4), kpts)):
763 raise ValueError(
764 'In explicit kpt list each row should have 4 elements')
766 clear_mp_keywords()
767 self.cell.kpoint_list = kpts
769 # Case 4: list or tuple of MP samples e.g. [3, 3, 2]
770 elif isinstance(kpts, (tuple, list)) and isinstance(kpts[0], int):
771 if len(kpts) != 3:
772 raise ValueError('Monkhorst-pack grid should have 3 values')
773 clear_mp_keywords()
774 self.cell.kpoint_mp_grid = '%d %d %d' % tuple(kpts)
776 # Case 5: str representation of Case 3 e.g. '3 3 2'
777 elif isinstance(kpts, str):
778 self.set_kpts([int(x) for x in kpts.split()])
780 # Case 6: dict of options e.g. {'size': (3, 3, 2), 'gamma': True}
781 # 'spacing' is allowed but transformed to 'density' to get mesh/offset
782 elif isinstance(kpts, dict):
783 kpts = kpts.copy()
785 if (kpts.get('spacing') is not None
786 and kpts.get('density') is not None):
787 raise ValueError(
788 'Cannot set kpts spacing and density simultaneously.')
789 else:
790 if kpts.get('spacing') is not None:
791 kpts = kpts.copy()
792 spacing = kpts.pop('spacing')
793 kpts['density'] = 1 / (2 * np.pi * spacing)
795 clear_mp_keywords()
796 size, offsets = kpts2sizeandoffsets(atoms=self.atoms, **kpts)
797 self.cell.kpoint_mp_grid = '%d %d %d' % tuple(size)
798 self.cell.kpoint_mp_offset = '%f %f %f' % tuple(offsets)
800 # Case 7: some other iterator. Try treating as a list:
801 elif hasattr(kpts, '__iter__'):
802 self.set_kpts(list(kpts))
804 # Otherwise, give up
805 else:
806 raise TypeError('Cannot interpret kpts of this type')
808 def todict(self, skip_default=True):
809 """Create dict with settings of .param and .cell"""
810 dct = {}
811 dct['param'] = self.param.get_attr_dict()
812 dct['cell'] = self.cell.get_attr_dict()
814 return dct
816 def check_state(self, atoms, tol=1e-15):
817 """Check for system changes since last calculation."""
818 return compare_atoms(self._old_atoms, atoms)
820 def _castep_find_last_record(self, castep_file):
821 """Checks wether a given castep file has a regular
822 ending message following the last banner message. If this
823 is the case, the line number of the last banner is message
824 is return, otherwise False.
826 returns (record_start, record_end, end_found, last_record_complete)
827 """
828 if isinstance(castep_file, str):
829 castep_file = paropen(castep_file, 'r')
830 file_opened = True
831 else:
832 file_opened = False
833 record_starts = []
834 while True:
835 line = castep_file.readline()
836 if (('Welcome' in line or 'Materials Studio' in line)
837 and 'CASTEP' in line):
838 record_starts = [castep_file.tell()] + record_starts
839 if not line:
840 break
842 if record_starts == []:
843 warnings.warn('Could not find CASTEP label in result file: %s.'
844 ' Are you sure this is a .castep file?' % castep_file)
845 return
847 # search for regular end of file
848 end_found = False
849 # start to search from record beginning from the back
850 # and see if
851 record_end = -1
852 for record_nr, record_start in enumerate(record_starts):
853 castep_file.seek(record_start)
854 while True:
855 line = castep_file.readline()
856 if not line:
857 break
858 if 'warn' in line.lower():
859 self._warnings.append(line)
861 if 'Finalisation time =' in line:
862 end_found = True
863 record_end = castep_file.tell()
864 break
866 if end_found:
867 break
869 if file_opened:
870 castep_file.close()
872 if end_found:
873 # record_nr == 0 corresponds to the last record here
874 if record_nr == 0:
875 return (record_start, record_end, True, True)
876 else:
877 return (record_start, record_end, True, False)
878 else:
879 return (0, record_end, False, False)
881 def read(self, castep_file=None):
882 """Read a castep file into the current instance."""
884 _close = True
886 if castep_file is None:
887 if self._castep_file:
888 castep_file = self._castep_file
889 out = paropen(castep_file, 'r')
890 else:
891 warnings.warn('No CASTEP file specified')
892 return
893 if not os.path.exists(castep_file):
894 warnings.warn('No CASTEP file found')
896 elif isinstance(castep_file, str):
897 out = paropen(castep_file, 'r')
899 else:
900 # in this case we assume that we have a fileobj already, but check
901 # for attributes in order to avoid extended EAFP blocks.
902 out = castep_file
904 # look before you leap...
905 attributes = ['name',
906 'seek',
907 'close',
908 'readline',
909 'tell']
911 for attr in attributes:
912 if not hasattr(out, attr):
913 raise TypeError(
914 '"castep_file" is neither str nor valid fileobj')
916 castep_file = out.name
917 _close = False
919 if self._seed is None:
920 self._seed = os.path.splitext(os.path.basename(castep_file))[0]
922 err_file = f'{self._seed}.0001.err'
923 if os.path.exists(err_file):
924 err_file = paropen(err_file)
925 self._error = err_file.read()
926 err_file.close()
927 # we return right-away because it might
928 # just be here from a previous run
929 # look for last result, if several CASTEP
930 # run are appended
932 record_start, record_end, end_found, _\
933 = self._castep_find_last_record(out)
934 if not end_found:
935 warnings.warn(
936 f'No regular end found in {castep_file} file. {self._error}')
937 if _close:
938 out.close()
939 return
940 # we return here, because the file has no a regular end
942 # now iterate over last CASTEP output in file to extract information
943 # could be generalized as well to extract trajectory from file
944 # holding several outputs
945 n_cell_const = 0
946 forces = []
948 # HOTFIX:
949 # we have to initialize the _stress variable as a zero array
950 # otherwise the calculator crashes upon pickling trajectories
951 # Alternative would be to raise a NotImplementedError() which
952 # is also kind of not true, since we can extract stresses if
953 # the user configures CASTEP to print them in the outfile
954 # stress = []
955 stress = np.zeros([3, 3])
956 hirsh_volrat = []
958 # Two flags to check whether spin-polarized or not, and whether
959 # Hirshfeld volumes are calculated
960 spin_polarized = False
961 calculate_hirshfeld = False
962 mulliken_analysis = False
963 hirshfeld_analysis = False
964 kpoints = None
966 positions_frac_list = []
967 mulliken_charges = []
968 spins = []
970 out.seek(record_start)
971 while True:
972 # TODO: add a switch if we have a geometry optimization: record
973 # atoms objects for intermediate steps.
974 try:
975 line = out.readline()
976 if not line or out.tell() > record_end:
977 break
978 elif 'Hirshfeld Analysis' in line:
979 hirshfeld_charges = []
981 hirshfeld_analysis = True
982 # skip the separating line
983 line = out.readline()
984 # this is the headline
985 line = out.readline()
987 if 'Charge' in line:
988 # skip the next separator line
989 line = out.readline()
990 while True:
991 line = out.readline()
992 fields = line.split()
993 if len(fields) == 1:
994 break
995 else:
996 hirshfeld_charges.append(float(fields[-1]))
997 elif 'stress calculation' in line:
998 if line.split()[-1].strip() == 'on':
999 self.param.calculate_stress = True
1000 elif 'basis set accuracy' in line:
1001 self.param.basis_precision = line.split()[-1]
1002 elif 'plane wave basis set cut-off' in line:
1003 # NB this is set as a private "result" attribute to avoid
1004 # conflict with input option basis_precision
1005 cutoff = float(line.split()[-2])
1006 self._cut_off_energy = cutoff
1007 if self.param.basis_precision.value is None:
1008 self.param.cut_off_energy = cutoff
1009 elif 'total energy / atom convergence tol.' in line:
1010 elec_energy_tol = float(line.split()[-2])
1011 self.param.elec_energy_tol = elec_energy_tol
1012 elif 'convergence tolerance window' in line:
1013 elec_convergence_win = int(line.split()[-2])
1014 self.param.elec_convergence_win = elec_convergence_win
1015 elif re.match(r'\sfinite basis set correction\s*:', line):
1016 finite_basis_corr = line.split()[-1]
1017 fbc_possibilities = {'none': 0,
1018 'manual': 1, 'automatic': 2}
1019 fbc = fbc_possibilities[finite_basis_corr]
1020 self.param.finite_basis_corr = fbc
1021 elif 'Treating system as non-metallic' in line:
1022 self.param.fix_occupancy = True
1023 elif 'max. number of SCF cycles:' in line:
1024 max_no_scf = float(line.split()[-1])
1025 self.param.max_scf_cycles = max_no_scf
1026 elif 'density-mixing scheme' in line:
1027 mixing_scheme = line.split()[-1]
1028 self.param.mixing_scheme = mixing_scheme
1029 elif 'dump wavefunctions every' in line:
1030 no_dump_cycles = float(line.split()[-3])
1031 self.param.num_dump_cycles = no_dump_cycles
1032 elif 'optimization strategy' in line:
1033 lspl = line.split(":")
1034 if lspl[0].strip() != 'optimization strategy':
1035 # This can happen in iprint: 3 calculations
1036 continue
1037 if 'memory' in line:
1038 self.param.opt_strategy = 'Memory'
1039 if 'speed' in line:
1040 self.param.opt_strategy = 'Speed'
1041 elif 'calculation limited to maximum' in line:
1042 calc_limit = float(line.split()[-2])
1043 self.param.run_time = calc_limit
1044 elif 'type of calculation' in line:
1045 lspl = line.split(":")
1046 if lspl[0].strip() != 'type of calculation':
1047 # This can happen in iprint: 3 calculations
1048 continue
1049 calc_type = lspl[-1]
1050 calc_type = re.sub(r'\s+', ' ', calc_type)
1051 calc_type = calc_type.strip()
1052 if calc_type != 'single point energy':
1053 calc_type_possibilities = {
1054 'geometry optimization': 'GeometryOptimization',
1055 'band structure': 'BandStructure',
1056 'molecular dynamics': 'MolecularDynamics',
1057 'optical properties': 'Optics',
1058 'phonon calculation': 'Phonon',
1059 'E-field calculation': 'Efield',
1060 'Phonon followed by E-field': 'Phonon+Efield',
1061 'transition state search': 'TransitionStateSearch',
1062 'Magnetic Resonance': 'MagRes',
1063 'Core level spectra': 'Elnes',
1064 'Electronic Spectroscopy': 'ElectronicSpectroscopy'
1065 }
1066 ctype = calc_type_possibilities[calc_type]
1067 self.param.task = ctype
1068 elif 'using functional' in line:
1069 used_functional = line.split(":")[-1]
1070 used_functional = re.sub(r'\s+', ' ', used_functional)
1071 used_functional = used_functional.strip()
1072 if used_functional != 'Local Density Approximation':
1073 used_functional_possibilities = {
1074 'Perdew Wang (1991)': 'PW91',
1075 'Perdew Burke Ernzerhof': 'PBE',
1076 'revised Perdew Burke Ernzerhof': 'RPBE',
1077 'PBE with Wu-Cohen exchange': 'WC',
1078 'PBE for solids (2008)': 'PBESOL',
1079 'Hartree-Fock': 'HF',
1080 'Hartree-Fock +': 'HF-LDA',
1081 'Screened Hartree-Fock': 'sX',
1082 'Screened Hartree-Fock + ': 'sX-LDA',
1083 'hybrid PBE0': 'PBE0',
1084 'hybrid B3LYP': 'B3LYP',
1085 'hybrid HSE03': 'HSE03',
1086 'hybrid HSE06': 'HSE06'
1087 }
1088 used_func = used_functional_possibilities[
1089 used_functional]
1090 self.param.xc_functional = used_func
1091 elif 'output verbosity' in line:
1092 iprint = int(line.split()[-1][1])
1093 if int(iprint) != 1:
1094 self.param.iprint = iprint
1095 elif 'treating system as spin-polarized' in line:
1096 spin_polarized = True
1097 self.param.spin_polarized = spin_polarized
1098 elif 'treating system as non-spin-polarized' in line:
1099 spin_polarized = False
1100 elif 'Number of kpoints used' in line:
1101 kpoints = int(line.split('=')[-1].strip())
1102 elif 'Unit Cell' in line:
1103 lattice_real = []
1104 lattice_reci = []
1105 while True:
1106 line = out.readline()
1107 fields = line.split()
1108 if len(fields) == 6:
1109 break
1110 for i in range(3):
1111 lattice_real.append([float(f) for f in fields[0:3]])
1112 lattice_reci.append([float(f) for f in fields[3:7]])
1113 line = out.readline()
1114 fields = line.split()
1115 elif 'Cell Contents' in line:
1116 while True:
1117 line = out.readline()
1118 if 'Total number of ions in cell' in line:
1119 n_atoms = int(line.split()[7])
1120 if 'Total number of species in cell' in line:
1121 int(line.split()[7])
1122 fields = line.split()
1123 if len(fields) == 0:
1124 break
1125 elif 'Fractional coordinates of atoms' in line:
1126 species = []
1127 custom_species = None # A CASTEP special thing
1128 positions_frac = []
1129 # positions_cart = []
1130 while True:
1131 line = out.readline()
1132 fields = line.split()
1133 if len(fields) == 7:
1134 break
1135 for n in range(n_atoms):
1136 spec_custom = fields[1].split(':', 1)
1137 elem = spec_custom[0]
1138 if len(spec_custom) > 1 and custom_species is None:
1139 # Add it to the custom info!
1140 custom_species = list(species)
1141 species.append(elem)
1142 if custom_species is not None:
1143 custom_species.append(fields[1])
1144 positions_frac.append([float(s) for s in fields[3:6]])
1145 line = out.readline()
1146 fields = line.split()
1147 positions_frac_list.append(positions_frac)
1148 elif 'Files used for pseudopotentials' in line:
1149 while True:
1150 line = out.readline()
1151 if 'Pseudopotential generated on-the-fly' in line:
1152 continue
1153 fields = line.split()
1154 if (len(fields) >= 2):
1155 elem, pp_file = fields
1156 self.cell.species_pot = (elem, pp_file)
1157 else:
1158 break
1159 elif 'k-Points For BZ Sampling' in line:
1160 # TODO: generalize for non-Monkhorst Pack case
1161 # (i.e. kpoint lists) -
1162 # kpoints_offset cannot be read this way and
1163 # is hence always set to None
1164 while True:
1165 line = out.readline()
1166 if not line.strip():
1167 break
1168 if 'MP grid size for SCF calculation' in line:
1169 # kpoints = ' '.join(line.split()[-3:])
1170 # self.kpoints_mp_grid = kpoints
1171 # self.kpoints_mp_offset = '0. 0. 0.'
1172 # not set here anymore because otherwise
1173 # two calculator objects go out of sync
1174 # after each calculation triggering unnecessary
1175 # recalculation
1176 break
1177 elif 'Number of cell constraints' in line:
1178 n_cell_const = int(line.split()[4])
1179 elif 'Final energy' in line:
1180 self._energy_total = float(line.split()[-2])
1181 elif 'Final free energy' in line:
1182 self._energy_free = float(line.split()[-2])
1183 elif 'NB est. 0K energy' in line:
1184 self._energy_0K = float(line.split()[-2])
1185 # check if we had a finite basis set correction
1186 elif 'Total energy corrected for finite basis set' in line:
1187 self._energy_total_corr = float(line.split()[-2])
1189 # Add support for dispersion correction
1190 # filtering due to SEDC is done in get_potential_energy
1191 elif 'Dispersion corrected final energy' in line:
1192 self._dispcorr_energy_total = float(line.split()[-2])
1193 elif 'Dispersion corrected final free energy' in line:
1194 self._dispcorr_energy_free = float(line.split()[-2])
1195 elif 'dispersion corrected est. 0K energy' in line:
1196 self._dispcorr_energy_0K = float(line.split()[-2])
1198 # remember to remove constraint labels in force components
1199 # (lacking a space behind the actual floating point number in
1200 # the CASTEP output)
1201 elif '******************** Forces *********************'\
1202 in line or\
1203 '************** Symmetrised Forces ***************'\
1204 in line or\
1205 '************** Constrained Symmetrised Forces ****'\
1206 '**********'\
1207 in line or\
1208 '******************** Constrained Forces **********'\
1209 '**********'\
1210 in line or\
1211 '******************* Unconstrained Forces *********'\
1212 '**********'\
1213 in line:
1214 fix = []
1215 fix_cart = []
1216 forces = []
1217 while True:
1218 line = out.readline()
1219 fields = line.split()
1220 if len(fields) == 7:
1221 break
1222 for n in range(n_atoms):
1223 consd = np.array([0, 0, 0])
1224 fxyz = [0, 0, 0]
1225 for (i, force_component) in enumerate(fields[-4:-1]):
1226 if force_component.count("(cons'd)") > 0:
1227 consd[i] = 1
1228 fxyz[i] = float(force_component.replace(
1229 "(cons'd)", ''))
1230 if consd.all():
1231 fix.append(n)
1232 elif consd.any():
1233 fix_cart.append(FixCartesian(n, consd))
1234 forces.append(fxyz)
1235 line = out.readline()
1236 fields = line.split()
1238 # add support for Hirshfeld analysis
1239 elif 'Hirshfeld / free atomic volume :' in line:
1240 # if we are here, then params must be able to cope with
1241 # Hirshfeld flag (if castep_keywords.py matches employed
1242 # castep version)
1243 calculate_hirshfeld = True
1244 hirsh_volrat = []
1245 while True:
1246 line = out.readline()
1247 fields = line.split()
1248 if len(fields) == 1:
1249 break
1250 for n in range(n_atoms):
1251 hirsh_atom = float(fields[0])
1252 hirsh_volrat.append(hirsh_atom)
1253 while True:
1254 line = out.readline()
1255 if 'Hirshfeld / free atomic volume :' in line or\
1256 'Hirshfeld Analysis' in line:
1257 break
1258 line = out.readline()
1259 fields = line.split()
1261 elif '***************** Stress Tensor *****************'\
1262 in line or\
1263 '*********** Symmetrised Stress Tensor ***********'\
1264 in line:
1265 stress = []
1266 while True:
1267 line = out.readline()
1268 fields = line.split()
1269 if len(fields) == 6:
1270 break
1271 for n in range(3):
1272 stress.append([float(s) for s in fields[2:5]])
1273 line = out.readline()
1274 fields = line.split()
1275 line = out.readline()
1276 if "Pressure:" in line:
1277 self._pressure = float(line.split()[-2]) * units.GPa
1278 elif ('BFGS: starting iteration' in line
1279 or 'BFGS: improving iteration' in line):
1280 if n_cell_const < 6:
1281 lattice_real = []
1282 lattice_reci = []
1283 # backup previous configuration first:
1284 # for highly symmetric systems (where essentially only the
1285 # stress is optimized, but the atomic positions) positions
1286 # are only printed once.
1287 if species:
1288 prev_species = deepcopy(species)
1289 if positions_frac:
1290 prev_positions_frac = deepcopy(positions_frac)
1291 species = []
1292 positions_frac = []
1293 forces = []
1295 # HOTFIX:
1296 # Same reason for the stress initialization as before
1297 # stress = []
1298 stress = np.zeros([3, 3])
1300 # extract info from the Mulliken analysis
1301 elif 'Atomic Populations' in line:
1302 # sometimes this appears twice in a castep file
1304 mulliken_analysis = True
1305 # skip the separating line
1306 line = out.readline()
1307 # this is the headline
1308 line = out.readline()
1310 if 'Charge' in line:
1311 # skip the next separator line
1312 line = out.readline()
1313 while True:
1314 line = out.readline()
1315 fields = line.split()
1316 if len(fields) == 1:
1317 break
1319 # the check for len==7 is due to CASTEP 18
1320 # outformat changes
1321 if spin_polarized:
1322 if len(fields) != 7:
1323 spins.append(float(fields[-1]))
1324 mulliken_charges.append(float(fields[-2]))
1325 else:
1326 mulliken_charges.append(float(fields[-1]))
1328 # There is actually no good reason to get out of the loop
1329 # already at this point... or do I miss something?
1330 # elif 'BFGS: Final Configuration:' in line:
1331 # break
1332 elif 'warn' in line.lower():
1333 self._warnings.append(line)
1335 # fetch some last info
1336 elif 'Total time' in line:
1337 pattern = r'.*=\s*([\d\.]+) s'
1338 self._total_time = float(re.search(pattern, line).group(1))
1340 elif 'Peak Memory Use' in line:
1341 pattern = r'.*=\s*([\d]+) kB'
1342 self._peak_memory = int(re.search(pattern, line).group(1))
1344 except Exception as exception:
1345 sys.stderr.write(line + '|-> line triggered exception: '
1346 + str(exception))
1347 raise
1349 if _close:
1350 out.close()
1352 # in highly summetric crystals, positions and symmetry are only printed
1353 # upon init, hence we here restore these original values
1354 if not positions_frac:
1355 positions_frac = prev_positions_frac
1356 if not species:
1357 species = prev_species
1359 if not spin_polarized:
1360 # set to zero spin if non-spin polarized calculation
1361 spins = np.zeros(len(positions_frac))
1362 elif len(spins) != len(positions_frac):
1363 warnings.warn('Spins could not be read for the atoms despite'
1364 ' spin-polarized calculation; spins will be ignored')
1365 spins = np.zeros(len(positions_frac))
1367 positions_frac_atoms = np.array(positions_frac)
1368 forces_atoms = np.array(forces)
1369 spins_atoms = np.array(spins)
1371 if mulliken_analysis:
1372 if len(mulliken_charges) != len(positions_frac):
1373 warnings.warn(
1374 'Mulliken charges could not be read for the atoms;'
1375 ' charges will be ignored')
1376 mulliken_charges = np.zeros(len(positions_frac))
1377 mulliken_charges_atoms = np.array(mulliken_charges)
1378 else:
1379 mulliken_charges_atoms = np.zeros(len(positions_frac))
1381 if hirshfeld_analysis:
1382 hirshfeld_charges_atoms = np.array(hirshfeld_charges)
1383 else:
1384 hirshfeld_charges_atoms = None
1386 if calculate_hirshfeld:
1387 hirsh_atoms = np.array(hirsh_volrat)
1388 else:
1389 hirsh_atoms = np.zeros_like(spins)
1391 if self.atoms and not self._set_atoms:
1392 # compensate for internal reordering of atoms by CASTEP
1393 # using the fact that the order is kept within each species
1395 # positions_frac_ase = self.atoms.get_scaled_positions(wrap=False)
1396 atoms_assigned = [False] * len(self.atoms)
1398 # positions_frac_castep_init = np.array(positions_frac_list[0])
1399 positions_frac_castep = np.array(positions_frac_list[-1])
1401 # species_castep = list(species)
1402 forces_castep = np.array(forces)
1403 hirsh_castep = np.array(hirsh_volrat)
1404 spins_castep = np.array(spins)
1405 mulliken_charges_castep = np.array(mulliken_charges_atoms)
1407 # go through the atoms position list and replace
1408 # with the corresponding one from the
1409 # castep file corresponding atomic number
1410 for iase in range(n_atoms):
1411 for icastep in range(n_atoms):
1412 if (species[icastep] == self.atoms[iase].symbol
1413 and not atoms_assigned[icastep]):
1414 positions_frac_atoms[iase] = \
1415 positions_frac_castep[icastep]
1416 forces_atoms[iase] = np.array(forces_castep[icastep])
1417 if iprint > 1 and calculate_hirshfeld:
1418 hirsh_atoms[iase] = np.array(hirsh_castep[icastep])
1419 if spin_polarized:
1420 # reordering not necessary in case all spins == 0
1421 spins_atoms[iase] = np.array(spins_castep[icastep])
1422 mulliken_charges_atoms[iase] = np.array(
1423 mulliken_charges_castep[icastep])
1424 atoms_assigned[icastep] = True
1425 break
1427 if not all(atoms_assigned):
1428 not_assigned = [i for (i, assigned)
1429 in zip(range(len(atoms_assigned)),
1430 atoms_assigned) if not assigned]
1431 warnings.warn(
1432 '%s atoms not assigned. '
1433 ' DEBUGINFO: The following atoms where not assigned: %s' %
1434 (atoms_assigned.count(False), not_assigned))
1435 else:
1436 self.atoms.set_scaled_positions(positions_frac_atoms)
1438 else:
1439 # If no atoms, object has been previously defined
1440 # we define it here and set the Castep() instance as calculator.
1441 # This covers the case that we simply want to open a .castep file.
1443 # The next time around we will have an atoms object, since
1444 # set_calculator also set atoms in the calculator.
1445 if self.atoms:
1446 constraints = self.atoms.constraints
1447 else:
1448 constraints = []
1449 atoms = ase.atoms.Atoms(species,
1450 cell=lattice_real,
1451 constraint=constraints,
1452 pbc=True,
1453 scaled_positions=positions_frac,
1454 )
1455 if custom_species is not None:
1456 atoms.new_array('castep_custom_species',
1457 np.array(custom_species))
1459 if self.param.spin_polarized:
1460 # only set magnetic moments if this was a spin polarized
1461 # calculation
1462 # this one fails as is
1463 atoms.set_initial_magnetic_moments(magmoms=spins_atoms)
1465 if mulliken_analysis:
1466 atoms.set_initial_charges(charges=mulliken_charges_atoms)
1467 atoms.calc = self
1469 self._kpoints = kpoints
1470 self._forces = forces_atoms
1471 # stress in .castep file is given in GPa:
1472 self._stress = np.array(stress) * units.GPa
1473 self._hirsh_volrat = hirsh_atoms
1474 self._spins = spins_atoms
1475 self._mulliken_charges = mulliken_charges_atoms
1476 self._hirshfeld_charges = hirshfeld_charges_atoms
1478 if self._warnings:
1479 warnings.warn(f'WARNING: {castep_file} contains warnings')
1480 for warning in self._warnings:
1481 warnings.warn(warning)
1482 # reset
1483 self._warnings = []
1485 # Read in eigenvalues from bands file
1486 bands_file = castep_file[:-7] + '.bands'
1487 if (self.param.task.value is not None
1488 and self.param.task.value.lower() == 'bandstructure'):
1489 self._band_structure = self.band_structure(bandfile=bands_file)
1490 else:
1491 try:
1492 (self._ibz_kpts,
1493 self._ibz_weights,
1494 self._eigenvalues,
1495 self._efermi) = read_bands(filename=bands_file)
1496 except FileNotFoundError:
1497 warnings.warn('Could not load .bands file, eigenvalues and '
1498 'Fermi energy are unknown')
1500 def get_hirsh_volrat(self):
1501 """
1502 Return the Hirshfeld volumes.
1503 """
1504 return self._hirsh_volrat
1506 def get_spins(self):
1507 """
1508 Return the spins from a plane-wave Mulliken analysis.
1509 """
1510 return self._spins
1512 def get_mulliken_charges(self):
1513 """
1514 Return the charges from a plane-wave Mulliken analysis.
1515 """
1516 return self._mulliken_charges
1518 def get_hirshfeld_charges(self):
1519 """
1520 Return the charges from a Hirshfeld analysis.
1521 """
1522 return self._hirshfeld_charges
1524 def get_total_time(self):
1525 """
1526 Return the total runtime
1527 """
1528 return self._total_time
1530 def get_peak_memory(self):
1531 """
1532 Return the peak memory usage
1533 """
1534 return self._peak_memory
1536 def set_label(self, label):
1537 """The label is part of each seed, which in turn is a prefix
1538 in each CASTEP related file.
1539 """
1540 # we may think about changing this in future to set `self._directory`
1541 # and `self._label`, as one would expect
1542 self._label = label
1544 def set_pspot(self, pspot, elems=None,
1545 notelems=None,
1546 clear=True,
1547 suffix='usp'):
1548 """Quickly set all pseudo-potentials: Usually CASTEP psp are named
1549 like <Elem>_<pspot>.<suffix> so this function function only expects
1550 the <LibraryName>. It then clears any previous pseudopotential
1551 settings apply the one with <LibraryName> for each element in the
1552 atoms object. The optional elems and notelems arguments can be used
1553 to exclusively assign to some species, or to exclude with notelemens.
1555 Parameters ::
1557 - elems (None) : set only these elements
1558 - notelems (None): do not set the elements
1559 - clear (True): clear previous settings
1560 - suffix (usp): PP file suffix
1561 """
1562 if self._find_pspots:
1563 if self._pedantic:
1564 warnings.warn('Warning: <_find_pspots> = True. '
1565 'Do you really want to use `set_pspots()`? '
1566 'This does not check whether the PP files exist. '
1567 'You may rather want to use `find_pspots()` with '
1568 'the same <pspot>.')
1570 if clear and not elems and not notelems:
1571 self.cell.species_pot.clear()
1572 for elem in set(self.atoms.get_chemical_symbols()):
1573 if elems is not None and elem not in elems:
1574 continue
1575 if notelems is not None and elem in notelems:
1576 continue
1577 self.cell.species_pot = (elem, f'{elem}_{pspot}.{suffix}')
1579 def find_pspots(self, pspot='.+', elems=None,
1580 notelems=None, clear=True, suffix='(usp|UPF|recpot)'):
1581 r"""Quickly find and set all pseudo-potentials by searching in
1582 castep_pp_path:
1584 This one is more flexible than set_pspots, and also checks if the files
1585 are actually available from the castep_pp_path.
1587 Essentially, the function parses the filenames in <castep_pp_path> and
1588 does a regex matching. The respective pattern is:
1590 r"^(<elem>|<elem.upper()>|elem.lower()>(_|-)<pspot>\.<suffix>$"
1592 In most cases, it will be sufficient to not specify anything, if you
1593 use standard CASTEP USPPs with only one file per element in the
1594 <castep_pp_path>.
1596 The function raises a `RuntimeError` if there is some ambiguity
1597 (multiple files per element).
1599 Parameters ::
1601 - pspots ('.+') : as defined above, will be a wildcard if not
1602 specified.
1603 - elems (None) : set only these elements
1604 - notelems (None): do not set the elements
1605 - clear (True): clear previous settings
1606 - suffix (usp|UPF|recpot): PP file suffix
1607 """
1608 if clear and not elems and not notelems:
1609 self.cell.species_pot.clear()
1611 if not os.path.isdir(self._castep_pp_path):
1612 if self._pedantic:
1613 warnings.warn(
1614 'Cannot search directory: {} Folder does not exist'
1615 .format(self._castep_pp_path))
1616 return
1618 # translate the bash wildcard syntax to regex
1619 if pspot == '*':
1620 pspot = '.*'
1621 if suffix == '*':
1622 suffix = '.*'
1623 if pspot == '*':
1624 pspot = '.*'
1626 # GBRV USPPs have a strnage naming schme
1627 pattern = r'^({elem}|{elem_upper}|{elem_lower})(_|-){pspot}\.{suffix}$'
1629 for elem in set(self.atoms.get_chemical_symbols()):
1630 if elems is not None and elem not in elems:
1631 continue
1632 if notelems is not None and elem in notelems:
1633 continue
1634 p = pattern.format(elem=elem,
1635 elem_upper=elem.upper(),
1636 elem_lower=elem.lower(),
1637 pspot=pspot,
1638 suffix=suffix)
1639 pps = []
1640 for f in os.listdir(self._castep_pp_path):
1641 if re.match(p, f):
1642 pps.append(f)
1643 if not pps:
1644 if self._pedantic:
1645 warnings.warn('Pseudopotential for species {} not found!'
1646 .format(elem))
1647 elif not len(pps) == 1:
1648 raise RuntimeError(
1649 'Pseudopotential for species ''{} not unique!\n'
1650 .format(elem)
1651 + 'Found the following files in {}\n'
1652 .format(self._castep_pp_path)
1653 + '\n'.join([f' {pp}' for pp in pps]) +
1654 '\nConsider a stricter search pattern in `find_pspots()`.')
1655 else:
1656 self.cell.species_pot = (elem, pps[0])
1658 @property
1659 def name(self):
1660 """Return the name of the calculator (string). """
1661 return self.__name__
1663 def get_property(self, name, atoms=None, allow_calculation=True):
1664 # High-level getter for compliance with the database module...
1665 # in principle this would not be necessary any longer if we properly
1666 # based this class on `Calculator`
1667 if name == 'forces':
1668 return self.get_forces(atoms)
1669 elif name == 'energy':
1670 return self.get_potential_energy(atoms)
1671 elif name == 'stress':
1672 return self.get_stress(atoms)
1673 elif name == 'charges':
1674 return self.get_charges(atoms)
1675 else:
1676 raise PropertyNotImplementedError
1678 @_self_getter
1679 def get_forces(self, atoms):
1680 """Run CASTEP calculation if needed and return forces."""
1681 self.update(atoms)
1682 return np.array(self._forces)
1684 @_self_getter
1685 def get_total_energy(self, atoms):
1686 """Run CASTEP calculation if needed and return total energy."""
1687 self.update(atoms)
1688 return self._energy_total
1690 @_self_getter
1691 def get_total_energy_corrected(self, atoms):
1692 """Run CASTEP calculation if needed and return total energy."""
1693 self.update(atoms)
1694 return self._energy_total_corr
1696 @_self_getter
1697 def get_free_energy(self, atoms):
1698 """Run CASTEP calculation if needed and return free energy.
1699 Only defined with smearing."""
1700 self.update(atoms)
1701 return self._energy_free
1703 @_self_getter
1704 def get_0K_energy(self, atoms):
1705 """Run CASTEP calculation if needed and return 0K energy.
1706 Only defined with smearing."""
1707 self.update(atoms)
1708 return self._energy_0K
1710 @_self_getter
1711 def get_potential_energy(self, atoms, force_consistent=False):
1712 # here for compatibility with ase/calculators/general.py
1713 # but accessing only _name variables
1714 """Return the total potential energy."""
1715 self.update(atoms)
1716 if force_consistent:
1717 # Assumption: If no dispersion correction is applied, then the
1718 # respective value will default to None as initialized.
1719 if self._dispcorr_energy_free is not None:
1720 return self._dispcorr_energy_free
1721 else:
1722 return self._energy_free
1723 else:
1724 if self._energy_0K is not None:
1725 if self._dispcorr_energy_0K is not None:
1726 return self._dispcorr_energy_0K
1727 else:
1728 return self._energy_0K
1729 else:
1730 if self._dispcorr_energy_total is not None:
1731 return self._dispcorr_energy_total
1732 else:
1733 if self._energy_total_corr is not None:
1734 return self._energy_total_corr
1735 else:
1736 return self._energy_total
1738 @_self_getter
1739 def get_stress(self, atoms):
1740 """Return the stress."""
1741 self.update(atoms)
1742 # modification: we return the Voigt form directly to get rid of the
1743 # annoying user warnings
1744 stress = np.array(
1745 [self._stress[0, 0], self._stress[1, 1], self._stress[2, 2],
1746 self._stress[1, 2], self._stress[0, 2], self._stress[0, 1]])
1747 # return self._stress
1748 return stress
1750 @_self_getter
1751 def get_pressure(self, atoms):
1752 """Return the pressure."""
1753 self.update(atoms)
1754 return self._pressure
1756 @_self_getter
1757 def get_unit_cell(self, atoms):
1758 """Return the unit cell."""
1759 self.update(atoms)
1760 return self._unit_cell
1762 @_self_getter
1763 def get_kpoints(self, atoms):
1764 """Return the kpoints."""
1765 self.update(atoms)
1766 return self._kpoints
1768 @_self_getter
1769 def get_number_cell_constraints(self, atoms):
1770 """Return the number of cell constraints."""
1771 self.update(atoms)
1772 return self._number_of_cell_constraints
1774 @_self_getter
1775 def get_charges(self, atoms):
1776 """Run CASTEP calculation if needed and return Mulliken charges."""
1777 self.update(atoms)
1778 return np.array(self._mulliken_charges)
1780 @_self_getter
1781 def get_magnetic_moments(self, atoms):
1782 """Run CASTEP calculation if needed and return Mulliken charges."""
1783 self.update(atoms)
1784 return np.array(self._spins)
1786 def set_atoms(self, atoms):
1787 """Sets the atoms for the calculator and vice versa."""
1788 atoms.pbc = [True, True, True]
1789 self.__dict__['atoms'] = atoms.copy()
1790 self.atoms._calc = self
1792 def update(self, atoms):
1793 """Checks if atoms object or calculator changed and
1794 runs calculation if so.
1795 """
1796 if self.calculation_required(atoms):
1797 self.calculate(atoms)
1799 def calculation_required(self, atoms, _=None):
1800 """Checks wether anything changed in the atoms object or CASTEP
1801 settings since the last calculation using this instance.
1802 """
1803 # SPR: what happens with the atoms parameter here? Why don't we use it?
1804 # from all that I can tell we need to compare against atoms instead of
1805 # self.atoms
1806 # if not self.atoms == self._old_atoms:
1807 if not atoms == self._old_atoms:
1808 return True
1809 if self._old_param is None or self._old_cell is None:
1810 return True
1811 if not self.param._options == self._old_param._options:
1812 return True
1813 if not self.cell._options == self._old_cell._options:
1814 return True
1815 return False
1817 def calculate(self, atoms):
1818 """Write all necessary input file and call CASTEP."""
1819 self.prepare_input_files(atoms, force_write=self._force_write)
1820 if not self._prepare_input_only:
1821 self.run()
1822 self.read()
1824 # we need to push the old state here!
1825 # although run() pushes it, read() may change the atoms object
1826 # again.
1827 # yet, the old state is supposed to be the one AFTER read()
1828 self.push_oldstate()
1830 def push_oldstate(self):
1831 """This function pushes the current state of the (CASTEP) Atoms object
1832 onto the previous state. Or in other words after calling this function,
1833 calculation_required will return False and enquiry functions just
1834 report the current value, e.g. get_forces(), get_potential_energy().
1835 """
1836 # make a snapshot of all current input
1837 # to be able to test if recalculation
1838 # is necessary
1839 self._old_atoms = self.atoms.copy()
1840 self._old_param = deepcopy(self.param)
1841 self._old_cell = deepcopy(self.cell)
1843 def initialize(self, *args, **kwargs):
1844 """Just an alias for prepar_input_files to comply with standard
1845 function names in ASE.
1846 """
1847 self.prepare_input_files(*args, **kwargs)
1849 def prepare_input_files(self, atoms=None, force_write=None):
1850 """Only writes the input .cell and .param files and return
1851 This can be useful if one quickly needs to prepare input files
1852 for a cluster where no python or ASE is available. One can than
1853 upload the file manually and read out the results using
1854 Castep().read().
1855 """
1857 if self.param.reuse.value is None:
1858 if self._pedantic:
1859 warnings.warn(
1860 'You have not set e.g. calc.param.reuse = True. '
1861 'Reusing a previous calculation may save CPU time! '
1862 'The interface will make sure by default, .check exists. '
1863 'file before adding this statement to the .param file.')
1864 if self.param.num_dump_cycles.value is None:
1865 if self._pedantic:
1866 warnings.warn(
1867 'You have not set e.g. calc.param.num_dump_cycles = 0. '
1868 'This can save you a lot of disk space. One only needs '
1869 '*wvfn* if electronic convergence is not achieved.')
1870 from ase.io.castep import write_param
1872 if atoms is None:
1873 atoms = self.atoms
1874 else:
1875 self.atoms = atoms
1877 if force_write is None:
1878 force_write = self._force_write
1880 # if we have new instance of the calculator,
1881 # move existing results out of the way, first
1882 if (os.path.isdir(self._directory)
1883 and self._calls == 0
1884 and self._rename_existing_dir):
1885 if os.listdir(self._directory) == []:
1886 os.rmdir(self._directory)
1887 else:
1888 # rename appending creation date of the directory
1889 ctime = time.localtime(os.lstat(self._directory).st_ctime)
1890 os.rename(self._directory, '%s.bak-%s' %
1891 (self._directory,
1892 time.strftime('%Y%m%d-%H%M%S', ctime)))
1894 # create work directory
1895 if not os.path.isdir(self._directory):
1896 os.makedirs(self._directory, 0o775)
1898 # we do this every time, not only upon first call
1899 # if self._calls == 0:
1900 self._fetch_pspots()
1902 # if _try_reuse is requested and this
1903 # is not the first run, we try to find
1904 # the .check file from the previous run
1905 # this is only necessary if _track_output
1906 # is set to true
1907 if self._try_reuse and self._calls > 0:
1908 if os.path.exists(self._abs_path(self._check_file)):
1909 self.param.reuse = self._check_file
1910 elif os.path.exists(self._abs_path(self._castep_bin_file)):
1911 self.param.reuse = self._castep_bin_file
1912 self._seed = self._build_castep_seed()
1913 self._check_file = f'{self._seed}.check'
1914 self._castep_bin_file = f'{self._seed}.castep_bin'
1915 self._castep_file = self._abs_path(f'{self._seed}.castep')
1917 # write out the input file
1918 self._write_cell(self._abs_path(f'{self._seed}.cell'),
1919 self.atoms, castep_cell=self.cell,
1920 force_write=force_write)
1922 if self._export_settings:
1923 interface_options = self._opt
1924 else:
1925 interface_options = None
1926 write_param(self._abs_path(f'{self._seed}.param'), self.param,
1927 check_checkfile=self._check_checkfile,
1928 force_write=force_write,
1929 interface_options=interface_options,)
1931 def _build_castep_seed(self):
1932 """Abstracts to construction of the final castep <seed>
1933 with and without _tracking_output.
1934 """
1935 if self._track_output:
1936 return '%s-%06d' % (self._label, self._calls)
1937 else:
1938 return f'{(self._label)}'
1940 def _abs_path(self, path):
1941 # Create an absolute path for a file to put in the working directory
1942 return os.path.join(self._directory, path)
1944 def run(self):
1945 """Simply call castep. If the first .err file
1946 contains text, this will be printed to the screen.
1947 """
1948 # change to target directory
1949 self._calls += 1
1951 # run castep itself
1952 stdout, stderr = shell_stdouterr('{} {}'.format(self._castep_command,
1953 self._seed),
1954 cwd=self._directory)
1955 if stdout:
1956 print(f'castep call stdout:\n{stdout}')
1957 if stderr:
1958 print(f'castep call stderr:\n{stderr}')
1960 # shouldn't it be called after read()???
1961 # self.push_oldstate()
1963 # check for non-empty error files
1964 err_file = self._abs_path(f'{self._seed}.0001.err')
1965 if os.path.exists(err_file):
1966 err_file = open(err_file)
1967 self._error = err_file.read()
1968 err_file.close()
1969 if self._error:
1970 raise RuntimeError(self._error)
1972 def __repr__(self):
1973 """Returns generic, fast to capture representation of
1974 CASTEP settings along with atoms object.
1975 """
1976 expr = ''
1977 expr += '-----------------Atoms--------------------\n'
1978 if self.atoms is not None:
1979 expr += str('%20s\n' % self.atoms)
1980 else:
1981 expr += 'None\n'
1983 expr += '-----------------Param keywords-----------\n'
1984 expr += str(self.param)
1985 expr += '-----------------Cell keywords------------\n'
1986 expr += str(self.cell)
1987 expr += '-----------------Internal keys------------\n'
1988 for key in self.internal_keys:
1989 expr += '%20s : %s\n' % (key, self._opt[key])
1991 return expr
1993 def __getattr__(self, attr):
1994 """___getattr___ gets overloaded to reroute the internal keys
1995 and to be able to easily store them in in the param so that
1996 they can be read in again in subsequent calls.
1997 """
1998 if attr in self.internal_keys:
1999 return self._opt[attr]
2000 if attr in ['__repr__', '__str__']:
2001 raise AttributeError
2002 elif attr not in self.__dict__:
2003 raise AttributeError(f'Attribute {attr} not found')
2004 else:
2005 return self.__dict__[attr]
2007 def __setattr__(self, attr, value):
2008 """We overload the settattr method to make value assignment
2009 as pythonic as possible. Internal values all start with _.
2010 Value assigment is case insensitive!
2011 """
2013 if attr.startswith('_'):
2014 # internal variables all start with _
2015 # let's check first if they are close but not identical
2016 # to one of the switches, that the user accesses directly
2017 similars = difflib.get_close_matches(attr, self.internal_keys,
2018 cutoff=0.9)
2019 if attr not in self.internal_keys and similars:
2020 warnings.warn(
2021 'Warning: You probably tried one of: '
2022 f'{similars} but typed {attr}')
2023 if attr in self.internal_keys:
2024 self._opt[attr] = value
2025 if attr == '_track_output':
2026 if value:
2027 self._try_reuse = True
2028 if self._pedantic:
2029 warnings.warn(
2030 'You switched _track_output on. This will '
2031 'consume a lot of disk-space. The interface '
2032 'also switched _try_reuse on, which will '
2033 'try to find the last check file. Set '
2034 '_try_reuse = False, if you need '
2035 'really separate calculations')
2036 elif '_try_reuse' in self._opt and self._try_reuse:
2037 self._try_reuse = False
2038 if self._pedantic:
2039 warnings.warn('_try_reuse is set to False, too')
2040 else:
2041 self.__dict__[attr] = value
2042 return
2043 elif attr in ['atoms', 'cell', 'param']:
2044 if value is not None:
2045 if attr == 'atoms' and not isinstance(value, ase.atoms.Atoms):
2046 raise TypeError(
2047 f'{value} is not an instance of ase.atoms.Atoms.')
2048 elif attr == 'cell' and not isinstance(value, CastepCell):
2049 raise TypeError(
2050 f'{value} is not an instance of CastepCell.')
2051 elif attr == 'param' and not isinstance(value, CastepParam):
2052 raise TypeError(
2053 f'{value} is not an instance of CastepParam.')
2054 # These 3 are accepted right-away, no matter what
2055 self.__dict__[attr] = value
2056 return
2057 elif attr in self.atoms_obj_keys:
2058 # keywords which clearly belong to the atoms object are
2059 # rerouted to go there
2060 self.atoms.__dict__[attr] = value
2061 return
2062 elif attr in self.atoms_keys:
2063 # CASTEP keywords that should go into the atoms object
2064 # itself are blocked
2065 warnings.warn('Ignoring setings of "%s", since this has to be set '
2066 'through the atoms object' % attr)
2067 return
2069 attr = attr.lower()
2070 if attr not in (list(self.cell._options.keys())
2071 + list(self.param._options.keys())):
2072 # what is left now should be meant to be a castep keyword
2073 # so we first check if it defined, and if not offer some error
2074 # correction
2075 if self._kw_tol == 0:
2076 similars = difflib.get_close_matches(
2077 attr,
2078 self.cell._options.keys() + self.param._options.keys())
2079 if similars:
2080 raise RuntimeError(
2081 f'Option "{attr}" not known! You mean "{similars[0]}"?')
2082 else:
2083 raise RuntimeError(f'Option "{attr}" is not known!')
2084 else:
2085 warnings.warn('Option "%s" is not known - please set any new'
2086 ' options directly in the .cell or .param '
2087 'objects' % attr)
2088 return
2090 # here we know it must go into one of the component param or cell
2091 # so we first determine which one
2092 if attr in self.param._options.keys():
2093 comp = 'param'
2094 elif attr in self.cell._options.keys():
2095 comp = 'cell'
2096 else:
2097 raise RuntimeError('Programming error: could not attach '
2098 'the keyword to an input file')
2100 self.__dict__[comp].__setattr__(attr, value)
2102 def merge_param(self, param, overwrite=True, ignore_internal_keys=False):
2103 """Parse a param file and merge it into the current parameters."""
2104 if isinstance(param, CastepParam):
2105 for key, option in param._options.items():
2106 if option.value is not None:
2107 self.param.__setattr__(key, option.value)
2108 return
2110 elif isinstance(param, str):
2111 param_file = open(param)
2112 _close = True
2114 else:
2115 # in this case we assume that we have a fileobj already, but check
2116 # for attributes in order to avoid extended EAFP blocks.
2117 param_file = param
2119 # look before you leap...
2120 attributes = ['name',
2121 'close'
2122 'readlines']
2124 for attr in attributes:
2125 if not hasattr(param_file, attr):
2126 raise TypeError('"param" is neither CastepParam nor str '
2127 'nor valid fileobj')
2129 param = param_file.name
2130 _close = False
2132 self, int_opts = read_param(fd=param_file, calc=self,
2133 get_interface_options=True)
2135 # Add the interface options
2136 for k, val in int_opts.items():
2137 if (k in self.internal_keys and not ignore_internal_keys):
2138 if val in _tf_table:
2139 val = _tf_table[val]
2140 self._opt[k] = val
2142 if _close:
2143 param_file.close()
2145 def dryrun_ok(self, dryrun_flag='-dryrun'):
2146 """Starts a CASTEP run with the -dryrun flag [default]
2147 in a temporary and check wether all variables are initialized
2148 correctly. This is recommended for every bigger simulation.
2149 """
2150 from ase.io.castep import write_param
2152 temp_dir = tempfile.mkdtemp()
2153 self._fetch_pspots(temp_dir)
2154 seed = 'dryrun'
2156 self._write_cell(os.path.join(temp_dir, f'{seed}.cell'),
2157 self.atoms, castep_cell=self.cell)
2158 # This part needs to be modified now that we rely on the new formats.py
2159 # interface
2160 if not os.path.isfile(os.path.join(temp_dir, f'{seed}.cell')):
2161 warnings.warn(f'{seed}.cell not written - aborting dryrun')
2162 return
2163 write_param(os.path.join(temp_dir, f'{seed}.param'), self.param, )
2165 stdout, stderr = shell_stdouterr(('{} {} {}'.format(
2166 self._castep_command,
2167 seed,
2168 dryrun_flag)),
2169 cwd=temp_dir)
2171 if stdout:
2172 print(stdout)
2173 if stderr:
2174 print(stderr)
2175 result_file = open(os.path.join(temp_dir, f'{seed}.castep'))
2177 txt = result_file.read()
2178 ok_string = r'.*DRYRUN finished.*No problems found with input files.*'
2179 match = re.match(ok_string, txt, re.DOTALL)
2181 m = re.search(r'Number of kpoints used =\s*([0-9]+)', txt)
2182 if m:
2183 self._kpoints = int(m.group(1))
2184 else:
2185 warnings.warn(
2186 'Couldn\'t fetch number of kpoints from dryrun CASTEP file')
2188 err_file = os.path.join(temp_dir, f'{seed}.0001.err')
2189 if match is None and os.path.exists(err_file):
2190 err_file = open(err_file)
2191 self._error = err_file.read()
2192 err_file.close()
2194 result_file.close()
2195 shutil.rmtree(temp_dir)
2197 # re.match return None is the string does not match
2198 return match is not None
2200 def _fetch_pspots(self, directory=None):
2201 """Put all specified pseudo-potentials into the working directory.
2202 """
2203 # should be a '==' right? Otherwise setting _castep_pp_path is not
2204 # honored.
2205 if (not cfg.get('PSPOT_DIR', None)
2206 and self._castep_pp_path == os.path.abspath('.')):
2207 # By default CASTEP consults the environment variable
2208 # PSPOT_DIR. If this contains a list of colon separated
2209 # directories it will check those directories for pseudo-
2210 # potential files if not in the current directory.
2211 # Thus if PSPOT_DIR is set there is nothing left to do.
2212 # If however PSPOT_DIR was been accidentally set
2213 # (e.g. with regards to a different program)
2214 # setting CASTEP_PP_PATH to an explicit value will
2215 # still be honored.
2216 return
2218 if directory is None:
2219 directory = self._directory
2220 if not os.path.isdir(self._castep_pp_path):
2221 warnings.warn(f'PSPs directory {self._castep_pp_path} not found')
2222 pspots = {}
2223 if self._find_pspots:
2224 self.find_pspots()
2225 if self.cell.species_pot.value is not None:
2226 for line in self.cell.species_pot.value.split('\n'):
2227 line = line.split()
2228 if line:
2229 pspots[line[0]] = line[1]
2230 for species in self.atoms.get_chemical_symbols():
2231 if not pspots or species not in pspots.keys():
2232 if self._build_missing_pspots:
2233 if self._pedantic:
2234 warnings.warn(
2235 'Warning: you have no PP specified for %s. '
2236 'CASTEP will now generate an on-the-fly '
2237 'potentials. '
2238 'For sake of numerical consistency and efficiency '
2239 'this is discouraged.' % species)
2240 else:
2241 raise RuntimeError(
2242 f'Warning: you have no PP specified for {species}.')
2243 if self.cell.species_pot.value:
2244 for (species, pspot) in pspots.items():
2245 orig_pspot_file = os.path.join(self._castep_pp_path, pspot)
2246 cp_pspot_file = os.path.join(directory, pspot)
2247 if (os.path.exists(orig_pspot_file)
2248 and not os.path.exists(cp_pspot_file)):
2249 if self._copy_pspots:
2250 shutil.copy(orig_pspot_file, directory)
2251 elif self._link_pspots:
2252 os.symlink(orig_pspot_file, cp_pspot_file)
2253 else:
2254 if self._pedantic:
2255 warnings.warn(ppwarning)
2258ppwarning = ('Warning: PP files have neither been '
2259 'linked nor copied to the working directory. Make '
2260 'sure to set the evironment variable PSPOT_DIR '
2261 'accordingly!')
2264def get_castep_version(castep_command):
2265 """This returns the version number as printed in the CASTEP banner.
2266 For newer CASTEP versions ( > 6.1) the --version command line option
2267 has been added; this will be attempted first.
2268 """
2269 import tempfile
2270 with tempfile.TemporaryDirectory() as temp_dir:
2271 return _get_castep_version(castep_command, temp_dir)
2274def _get_castep_version(castep_command, temp_dir):
2275 jname = 'dummy_jobname'
2276 stdout, stderr = '', ''
2277 fallback_version = 16. # CASTEP 16.0 and 16.1 report version wrongly
2278 try:
2279 stdout, stderr = subprocess.Popen(
2280 castep_command.split() + ['--version'],
2281 stderr=subprocess.PIPE,
2282 stdout=subprocess.PIPE, cwd=temp_dir,
2283 universal_newlines=True).communicate()
2284 if 'CASTEP version' not in stdout:
2285 stdout, stderr = subprocess.Popen(
2286 castep_command.split() + [jname],
2287 stderr=subprocess.PIPE,
2288 stdout=subprocess.PIPE, cwd=temp_dir,
2289 universal_newlines=True).communicate()
2290 except Exception: # XXX Which kind of exception?
2291 msg = ''
2292 msg += 'Could not determine the version of your CASTEP binary \n'
2293 msg += 'This usually means one of the following \n'
2294 msg += ' * you do not have CASTEP installed \n'
2295 msg += ' * you have not set the CASTEP_COMMAND to call it \n'
2296 msg += ' * you have provided a wrong CASTEP_COMMAND. \n'
2297 msg += ' Make sure it is in your PATH\n\n'
2298 msg += stdout
2299 msg += stderr
2300 raise CastepVersionError(msg)
2301 if 'CASTEP version' in stdout:
2302 output_txt = stdout.split('\n')
2303 version_re = re.compile(r'CASTEP version:\s*([0-9\.]*)')
2304 else:
2305 output = open(os.path.join(temp_dir, f'{jname}.castep'))
2306 output_txt = output.readlines()
2307 output.close()
2308 version_re = re.compile(r'(?<=CASTEP version )[0-9.]*')
2309 # shutil.rmtree(temp_dir)
2310 for line in output_txt:
2311 if 'CASTEP version' in line:
2312 try:
2313 return float(version_re.findall(line)[0])
2314 except ValueError:
2315 # Fallback for buggy --version on CASTEP 16.0, 16.1
2316 return fallback_version
2319def create_castep_keywords(castep_command, filename='castep_keywords.json',
2320 force_write=True, path='.', fetch_only=None):
2321 """This function allows to fetch all available keywords from stdout
2322 of an installed castep binary. It furthermore collects the documentation
2323 to harness the power of (ipython) inspection and type for some basic
2324 type checking of input. All information is stored in a JSON file that is
2325 not distributed by default to avoid breaking the license of CASTEP.
2326 """
2327 # Takes a while ...
2328 # Fetch all allowed parameters
2329 # fetch_only : only fetch that many parameters (for testsuite only)
2330 suffixes = ['cell', 'param']
2332 filepath = os.path.join(path, filename)
2334 if os.path.exists(filepath) and not force_write:
2335 warnings.warn('CASTEP Options Module file exists. '
2336 'You can overwrite it by calling '
2337 'python castep.py -f [CASTEP_COMMAND].')
2338 return False
2340 # Not saving directly to file her to prevent half-generated files
2341 # which will cause problems on future runs
2343 castep_version = get_castep_version(castep_command)
2345 help_all, _ = shell_stdouterr(f'{castep_command} -help all')
2347 # Filter out proper keywords
2348 try:
2349 # The old pattern does not math properly as in CASTEP as of v8.0 there
2350 # are some keywords for the semi-empircal dispersion correction (SEDC)
2351 # which also include numbers.
2352 if castep_version < 7.0:
2353 pattern = r'((?<=^ )[A-Z_]{2,}|(?<=^)[A-Z_]{2,})'
2354 else:
2355 pattern = r'((?<=^ )[A-Z_\d]{2,}|(?<=^)[A-Z_\d]{2,})'
2357 raw_options = re.findall(pattern, help_all, re.MULTILINE)
2358 except Exception:
2359 warnings.warn(f'Problem parsing: {help_all}')
2360 raise
2362 types = set()
2363 levels = set()
2365 processed_n = 0
2366 to_process = len(raw_options[:fetch_only])
2368 processed_options = {sf: {} for sf in suffixes}
2370 for o_i, option in enumerate(raw_options[:fetch_only]):
2371 doc, _ = shell_stdouterr(f'{castep_command} -help {option}')
2373 # Stand Back! I know regular expressions (http://xkcd.com/208/) :-)
2374 match = re.match(r'(?P<before_type>.*)Type: (?P<type>.+?)\s+'
2375 + r'Level: (?P<level>[^ ]+)\n\s*\n'
2376 + r'(?P<doc>.*?)(\n\s*\n|$)', doc, re.DOTALL)
2378 processed_n += 1
2380 if match is not None:
2381 match = match.groupdict()
2383 # JM: uncomment lines in following block to debug issues
2384 # with keyword assignment during extraction process from CASTEP
2385 suffix = None
2386 if re.findall(r'PARAMETERS keywords:\n\n\s?None found', doc):
2387 suffix = 'cell'
2388 if re.findall(r'CELL keywords:\n\n\s?None found', doc):
2389 suffix = 'param'
2390 if suffix is None:
2391 warnings.warn('%s -> not assigned to either'
2392 ' CELL or PARAMETERS keywords' % option)
2394 option = option.lower()
2395 mtyp = match.get('type', None)
2396 mlvl = match.get('level', None)
2397 mdoc = match.get('doc', None)
2399 if mtyp is None:
2400 warnings.warn(f'Found no type for {option}')
2401 continue
2402 if mlvl is None:
2403 warnings.warn(f'Found no level for {option}')
2404 continue
2405 if mdoc is None:
2406 warnings.warn(f'Found no doc string for {option}')
2407 continue
2409 types = types.union([mtyp])
2410 levels = levels.union([mlvl])
2412 processed_options[suffix][option] = {
2413 'keyword': option,
2414 'option_type': mtyp,
2415 'level': mlvl,
2416 'docstring': mdoc
2417 }
2419 processed_n += 1
2421 frac = (o_i + 1.0) / to_process
2422 sys.stdout.write('\rProcessed: [{}] {:>3.0f}%'.format(
2423 '#' * int(frac * 20) + ' '
2424 * (20 - int(frac * 20)),
2425 100 * frac))
2426 sys.stdout.flush()
2428 else:
2429 warnings.warn(f'create_castep_keywords: Could not process {option}')
2431 sys.stdout.write('\n')
2432 sys.stdout.flush()
2434 processed_options['types'] = list(types)
2435 processed_options['levels'] = list(levels)
2436 processed_options['castep_version'] = castep_version
2438 json.dump(processed_options, open(filepath, 'w'), indent=4)
2440 warnings.warn(f'CASTEP v{castep_version}, fetched {processed_n} keywords')
2441 return True
2444class CastepOption:
2445 """"A CASTEP option. It handles basic conversions from string to its value
2446 type."""
2448 default_convert_types = {
2449 'boolean (logical)': 'bool',
2450 'defined': 'bool',
2451 'string': 'str',
2452 'integer': 'int',
2453 'real': 'float',
2454 'integer vector': 'int_vector',
2455 'real vector': 'float_vector',
2456 'physical': 'float_physical',
2457 'block': 'block'
2458 }
2460 def __init__(self, keyword, level, option_type, value=None,
2461 docstring='No information available'):
2462 self.keyword = keyword
2463 self.level = level
2464 self.type = option_type
2465 self._value = value
2466 self.__doc__ = docstring
2468 @property
2469 def value(self):
2471 if self._value is not None:
2472 if self.type.lower() in ('integer vector', 'real vector',
2473 'physical'):
2474 return ' '.join(map(str, self._value))
2475 elif self.type.lower() in ('boolean (logical)', 'defined'):
2476 return str(self._value).upper()
2477 else:
2478 return str(self._value)
2480 @property
2481 def raw_value(self):
2482 # The value, not converted to a string
2483 return self._value
2485 @value.setter # type: ignore[attr-defined, no-redef]
2486 def value(self, val):
2488 if val is None:
2489 self.clear()
2490 return
2492 ctype = self.default_convert_types.get(self.type.lower(), 'str')
2493 typeparse = f'_parse_{ctype}'
2494 try:
2495 self._value = getattr(self, typeparse)(val)
2496 except ValueError:
2497 raise ConversionError(ctype, self.keyword, val)
2499 def clear(self):
2500 """Reset the value of the option to None again"""
2501 self._value = None
2503 @staticmethod
2504 def _parse_bool(value):
2505 try:
2506 value = _tf_table[str(value).strip().title()]
2507 except (KeyError, ValueError):
2508 raise ValueError()
2509 return value
2511 @staticmethod
2512 def _parse_str(value):
2513 value = str(value)
2514 return value
2516 @staticmethod
2517 def _parse_int(value):
2518 value = int(value)
2519 return value
2521 @staticmethod
2522 def _parse_float(value):
2523 value = float(value)
2524 return value
2526 @staticmethod
2527 def _parse_int_vector(value):
2528 # Accepts either a string or an actual list/numpy array of ints
2529 if isinstance(value, str):
2530 if ',' in value:
2531 value = value.replace(',', ' ')
2532 value = list(map(int, value.split()))
2534 value = np.array(value)
2536 if value.shape != (3,) or value.dtype != int:
2537 raise ValueError()
2539 return list(value)
2541 @staticmethod
2542 def _parse_float_vector(value):
2543 # Accepts either a string or an actual list/numpy array of floats
2544 if isinstance(value, str):
2545 if ',' in value:
2546 value = value.replace(',', ' ')
2547 value = list(map(float, value.split()))
2549 value = np.array(value) * 1.0
2551 if value.shape != (3,) or value.dtype != float:
2552 raise ValueError()
2554 return list(value)
2556 @staticmethod
2557 def _parse_float_physical(value):
2558 # If this is a string containing units, saves them
2559 if isinstance(value, str):
2560 value = value.split()
2562 try:
2563 l = len(value)
2564 except TypeError:
2565 l = 1
2566 value = [value]
2568 if l == 1:
2569 try:
2570 value = (float(value[0]), '')
2571 except (TypeError, ValueError):
2572 raise ValueError()
2573 elif l == 2:
2574 try:
2575 value = (float(value[0]), value[1])
2576 except (TypeError, ValueError, IndexError):
2577 raise ValueError()
2578 else:
2579 raise ValueError()
2581 return value
2583 @staticmethod
2584 def _parse_block(value):
2586 if isinstance(value, str):
2587 return value
2588 elif hasattr(value, '__getitem__'):
2589 return '\n'.join(value) # Arrays of lines
2590 else:
2591 raise ValueError()
2593 def __repr__(self):
2594 if self._value:
2595 expr = ('Option: {keyword}({type}, {level}):\n{_value}\n'
2596 ).format(**self.__dict__)
2597 else:
2598 expr = ('Option: {keyword}[unset]({type}, {level})'
2599 ).format(**self.__dict__)
2600 return expr
2602 def __eq__(self, other):
2603 if not isinstance(other, CastepOption):
2604 return False
2605 else:
2606 return self.__dict__ == other.__dict__
2609class CastepOptionDict:
2610 """A dictionary-like object to hold a set of options for .cell or .param
2611 files loaded from a dictionary, for the sake of validation.
2613 Replaces the old CastepCellDict and CastepParamDict that were defined in
2614 the castep_keywords.py file.
2615 """
2617 def __init__(self, options=None):
2618 object.__init__(self)
2619 self._options = {} # ComparableDict is not needed any more as
2620 # CastepOptions can be compared directly now
2621 for kw in options:
2622 opt = CastepOption(**options[kw])
2623 self._options[opt.keyword] = opt
2624 self.__dict__[opt.keyword] = opt
2627class CastepInputFile:
2629 """Master class for CastepParam and CastepCell to inherit from"""
2631 _keyword_conflicts: List[Set[str]] = []
2633 def __init__(self, options_dict=None, keyword_tolerance=1):
2634 object.__init__(self)
2636 if options_dict is None:
2637 options_dict = CastepOptionDict({})
2639 self._options = options_dict._options
2640 self.__dict__.update(self._options)
2641 # keyword_tolerance means how strict the checks on new attributes are
2642 # 0 = no new attributes allowed
2643 # 1 = new attributes allowed, warning given
2644 # 2 = new attributes allowed, silent
2645 self._perm = np.clip(keyword_tolerance, 0, 2)
2647 # Compile a dictionary for quick check of conflict sets
2648 self._conflict_dict = {
2649 kw: set(cset).difference({kw})
2650 for cset in self._keyword_conflicts for kw in cset}
2652 def __repr__(self):
2653 expr = ''
2654 is_default = True
2655 for key, option in sorted(self._options.items()):
2656 if option.value is not None:
2657 is_default = False
2658 expr += ('%20s : %s\n' % (key, option.value))
2659 if is_default:
2660 expr = 'Default\n'
2662 expr += f'Keyword tolerance: {self._perm}'
2663 return expr
2665 def __setattr__(self, attr, value):
2667 # Hidden attributes are treated normally
2668 if attr.startswith('_'):
2669 self.__dict__[attr] = value
2670 return
2672 if attr not in self._options.keys():
2674 if self._perm > 0:
2675 # Do we consider it a string or a block?
2676 is_str = isinstance(value, str)
2677 is_block = False
2678 if ((hasattr(value, '__getitem__') and not is_str)
2679 or (is_str and len(value.split('\n')) > 1)):
2680 is_block = True
2682 if self._perm == 0:
2683 similars = difflib.get_close_matches(attr,
2684 self._options.keys())
2685 if similars:
2686 raise RuntimeError(
2687 f'Option "{attr}" not known! You mean "{similars[0]}"?')
2688 else:
2689 raise RuntimeError(f'Option "{attr}" is not known!')
2690 elif self._perm == 1:
2691 warnings.warn(('Option "%s" is not known and will '
2692 'be added as a %s') % (attr,
2693 ('block' if is_block else
2694 'string')))
2695 attr = attr.lower()
2696 opt = CastepOption(keyword=attr, level='Unknown',
2697 option_type='block' if is_block else 'string')
2698 self._options[attr] = opt
2699 self.__dict__[attr] = opt
2700 else:
2701 attr = attr.lower()
2702 opt = self._options[attr]
2704 if not opt.type.lower() == 'block' and isinstance(value, str):
2705 value = value.replace(':', ' ')
2707 # If it is, use the appropriate parser, unless a custom one is defined
2708 attrparse = f'_parse_{attr.lower()}'
2710 # Check for any conflicts if the value is not None
2711 if value is not None:
2712 cset = self._conflict_dict.get(attr.lower(), {})
2713 for c in cset:
2714 if (c in self._options and self._options[c].value):
2715 warnings.warn(
2716 'option "{attr}" conflicts with "{conflict}" in '
2717 'calculator. Setting "{conflict}" to '
2718 'None.'.format(attr=attr, conflict=c))
2719 self._options[c].value = None
2721 if hasattr(self, attrparse):
2722 self._options[attr].value = self.__getattribute__(attrparse)(value)
2723 else:
2724 self._options[attr].value = value
2726 def __getattr__(self, name):
2727 if name[0] == '_' or self._perm == 0:
2728 raise AttributeError()
2730 if self._perm == 1:
2731 warnings.warn(f'Option {(name)} is not known, returning None')
2733 return CastepOption(keyword='none', level='Unknown',
2734 option_type='string', value=None)
2736 def get_attr_dict(self, raw=False, types=False):
2737 """Settings that go into .param file in a traditional dict"""
2739 attrdict = {k: o.raw_value if raw else o.value
2740 for k, o in self._options.items() if o.value is not None}
2742 if types:
2743 for key, val in attrdict.items():
2744 attrdict[key] = (val, self._options[key].type)
2746 return attrdict
2749class CastepParam(CastepInputFile):
2750 """CastepParam abstracts the settings that go into the .param file"""
2752 _keyword_conflicts = [{'cut_off_energy', 'basis_precision'}, ]
2754 def __init__(self, castep_keywords, keyword_tolerance=1):
2755 self._castep_version = castep_keywords.castep_version
2756 CastepInputFile.__init__(self, castep_keywords.CastepParamDict(),
2757 keyword_tolerance)
2759 @property
2760 def castep_version(self):
2761 return self._castep_version
2763 # .param specific parsers
2764 def _parse_reuse(self, value):
2765 if value is None:
2766 return None # Reset the value
2767 try:
2768 if self._options['continuation'].value:
2769 warnings.warn('Cannot set reuse if continuation is set, and '
2770 'vice versa. Set the other to None, if you want '
2771 'this setting.')
2772 return None
2773 except KeyError:
2774 pass
2775 return 'default' if (value is True) else str(value)
2777 def _parse_continuation(self, value):
2778 if value is None:
2779 return None # Reset the value
2780 try:
2781 if self._options['reuse'].value:
2782 warnings.warn('Cannot set reuse if continuation is set, and '
2783 'vice versa. Set the other to None, if you want '
2784 'this setting.')
2785 return None
2786 except KeyError:
2787 pass
2788 return 'default' if (value is True) else str(value)
2791class CastepCell(CastepInputFile):
2793 """CastepCell abstracts all setting that go into the .cell file"""
2795 _keyword_conflicts = [
2796 {'kpoint_mp_grid', 'kpoint_mp_spacing', 'kpoint_list',
2797 'kpoints_mp_grid', 'kpoints_mp_spacing', 'kpoints_list'},
2798 {'bs_kpoint_mp_grid',
2799 'bs_kpoint_mp_spacing',
2800 'bs_kpoint_list',
2801 'bs_kpoint_path',
2802 'bs_kpoints_mp_grid',
2803 'bs_kpoints_mp_spacing',
2804 'bs_kpoints_list',
2805 'bs_kpoints_path'},
2806 {'spectral_kpoint_mp_grid',
2807 'spectral_kpoint_mp_spacing',
2808 'spectral_kpoint_list',
2809 'spectral_kpoint_path',
2810 'spectral_kpoints_mp_grid',
2811 'spectral_kpoints_mp_spacing',
2812 'spectral_kpoints_list',
2813 'spectral_kpoints_path'},
2814 {'phonon_kpoint_mp_grid',
2815 'phonon_kpoint_mp_spacing',
2816 'phonon_kpoint_list',
2817 'phonon_kpoint_path',
2818 'phonon_kpoints_mp_grid',
2819 'phonon_kpoints_mp_spacing',
2820 'phonon_kpoints_list',
2821 'phonon_kpoints_path'},
2822 {'fine_phonon_kpoint_mp_grid',
2823 'fine_phonon_kpoint_mp_spacing',
2824 'fine_phonon_kpoint_list',
2825 'fine_phonon_kpoint_path'},
2826 {'magres_kpoint_mp_grid',
2827 'magres_kpoint_mp_spacing',
2828 'magres_kpoint_list',
2829 'magres_kpoint_path'},
2830 {'elnes_kpoint_mp_grid',
2831 'elnes_kpoint_mp_spacing',
2832 'elnes_kpoint_list',
2833 'elnes_kpoint_path'},
2834 {'optics_kpoint_mp_grid',
2835 'optics_kpoint_mp_spacing',
2836 'optics_kpoint_list',
2837 'optics_kpoint_path'},
2838 {'supercell_kpoint_mp_grid',
2839 'supercell_kpoint_mp_spacing',
2840 'supercell_kpoint_list',
2841 'supercell_kpoint_path'}, ]
2843 def __init__(self, castep_keywords, keyword_tolerance=1):
2844 self._castep_version = castep_keywords.castep_version
2845 CastepInputFile.__init__(self, castep_keywords.CastepCellDict(),
2846 keyword_tolerance)
2848 @property
2849 def castep_version(self):
2850 return self._castep_version
2852 # .cell specific parsers
2853 def _parse_species_pot(self, value):
2855 # Single tuple
2856 if isinstance(value, tuple) and len(value) == 2:
2857 value = [value]
2858 # List of tuples
2859 if hasattr(value, '__getitem__'):
2860 pspots = [tuple(map(str.strip, x)) for x in value]
2861 if not all(map(lambda x: len(x) == 2, value)):
2862 warnings.warn(
2863 'Please specify pseudopotentials in python as '
2864 'a tuple or a list of tuples formatted like: '
2865 '(species, file), e.g. ("O", "path-to/O_OTFG.usp") '
2866 'Anything else will be ignored')
2867 return None
2869 text_block = self._options['species_pot'].value
2871 text_block = text_block if text_block else ''
2872 # Remove any duplicates
2873 for pp in pspots:
2874 text_block = re.sub(fr'\n?\s*{pp[0]}\s+.*', '', text_block)
2875 if pp[1]:
2876 text_block += '\n%s %s' % pp
2878 return text_block
2880 def _parse_symmetry_ops(self, value):
2881 if not isinstance(value, tuple) \
2882 or not len(value) == 2 \
2883 or not value[0].shape[1:] == (3, 3) \
2884 or not value[1].shape[1:] == (3,) \
2885 or not value[0].shape[0] == value[1].shape[0]:
2886 warnings.warn('Invalid symmetry_ops block, skipping')
2887 return
2888 # Now on to print...
2889 text_block = ''
2890 for op_i, (op_rot, op_tranls) in enumerate(zip(*value)):
2891 text_block += '\n'.join([' '.join([str(x) for x in row])
2892 for row in op_rot])
2893 text_block += '\n'
2894 text_block += ' '.join([str(x) for x in op_tranls])
2895 text_block += '\n\n'
2897 return text_block
2899 def _parse_positions_abs_intermediate(self, value):
2900 return _parse_tss_block(value)
2902 def _parse_positions_abs_product(self, value):
2903 return _parse_tss_block(value)
2905 def _parse_positions_frac_intermediate(self, value):
2906 return _parse_tss_block(value, True)
2908 def _parse_positions_frac_product(self, value):
2909 return _parse_tss_block(value, True)
2912CastepKeywords = namedtuple('CastepKeywords',
2913 ['CastepParamDict', 'CastepCellDict',
2914 'types', 'levels', 'castep_version'])
2916# We keep this just for naming consistency with older versions
2919def make_cell_dict(data=None):
2921 data = data if data is not None else {}
2923 class CastepCellDict(CastepOptionDict):
2924 def __init__(self):
2925 CastepOptionDict.__init__(self, data)
2927 return CastepCellDict
2930def make_param_dict(data=None):
2932 data = data if data is not None else {}
2934 class CastepParamDict(CastepOptionDict):
2935 def __init__(self):
2936 CastepOptionDict.__init__(self, data)
2938 return CastepParamDict
2941class CastepVersionError(Exception):
2942 """No special behaviour, works to signal when Castep can not be found"""
2945class ConversionError(Exception):
2947 """Print customized error for options that are not converted correctly
2948 and point out that they are maybe not implemented, yet"""
2950 def __init__(self, key_type, attr, value):
2951 Exception.__init__(self)
2952 self.key_type = key_type
2953 self.value = value
2954 self.attr = attr
2956 def __str__(self):
2957 return f'Could not convert {self.attr} = {self.value} '\
2958 + 'to {self.key_type}\n' \
2959 + 'This means you either tried to set a value of the wrong\n'\
2960 + 'type or this keyword needs some special care. Please feel\n'\
2961 + 'to add it to the corresponding __setattr__ method and send\n'\
2962 + f'the patch to {(contact_email)}, so we can all benefit.'
2965def get_castep_pp_path(castep_pp_path=''):
2966 """Abstract the quest for a CASTEP PSP directory."""
2967 if castep_pp_path:
2968 return os.path.abspath(os.path.expanduser(castep_pp_path))
2969 elif 'PSPOT_DIR' in cfg:
2970 return cfg['PSPOT_DIR']
2971 elif 'CASTEP_PP_PATH' in cfg:
2972 return cfg['CASTEP_PP_PATH']
2973 else:
2974 return os.path.abspath('.')
2977def get_castep_command(castep_command=''):
2978 """Abstract the quest for a castep_command string."""
2979 if castep_command:
2980 return castep_command
2981 elif 'CASTEP_COMMAND' in cfg:
2982 return cfg['CASTEP_COMMAND']
2983 else:
2984 return 'castep'
2987def shell_stdouterr(raw_command, cwd=None):
2988 """Abstracts the standard call of the commandline, when
2989 we are only interested in the stdout and stderr
2990 """
2991 stdout, stderr = subprocess.Popen(raw_command,
2992 stdout=subprocess.PIPE,
2993 stderr=subprocess.PIPE,
2994 universal_newlines=True,
2995 shell=True, cwd=cwd).communicate()
2996 return stdout.strip(), stderr.strip()
2999def import_castep_keywords(castep_command='',
3000 filename='castep_keywords.json',
3001 path='.'):
3002 """Search for castep keywords JSON in multiple paths"""
3004 config_paths = ('~/.ase', '~/.config/ase')
3005 searchpaths = [path] + [os.path.expanduser(config_path)
3006 for config_path in config_paths]
3007 try:
3008 keywords_file = sum([glob.glob(os.path.join(sp, filename))
3009 for sp in searchpaths], [])[0]
3010 except IndexError:
3011 warnings.warn("""Generating CASTEP keywords JSON file... hang on.
3012 The CASTEP keywords JSON file contains abstractions for CASTEP input
3013 parameters (for both .cell and .param input files), including some
3014 format checks and descriptions. The latter are extracted from the
3015 internal online help facility of a CASTEP binary, thus allowing to
3016 easily keep the calculator synchronized with (different versions of)
3017 the CASTEP code. Consequently, avoiding licensing issues (CASTEP is
3018 distributed commercially by Biovia), we consider it wise not to
3019 provide the file in the first place.""")
3020 create_castep_keywords(get_castep_command(castep_command),
3021 filename=filename, path=path)
3022 keywords_file = Path(path).absolute() / filename
3024 warnings.warn(
3025 f'Stored castep keywords dictionary as {keywords_file}. '
3026 f'Copy it to {Path(config_paths[0]).expanduser() / filename} for '
3027 r'user installation.')
3029 # Now create the castep_keywords object proper
3030 with open(keywords_file) as fd:
3031 kwdata = json.load(fd)
3033 # This is a bit awkward, but it's necessary for backwards compatibility
3034 param_dict = make_param_dict(kwdata['param'])
3035 cell_dict = make_cell_dict(kwdata['cell'])
3037 castep_keywords = CastepKeywords(param_dict, cell_dict,
3038 kwdata['types'], kwdata['levels'],
3039 kwdata['castep_version'])
3041 return castep_keywords
3044if __name__ == '__main__':
3045 warnings.warn(
3046 'When called directly this calculator will fetch all available '
3047 'keywords from the binarys help function into a '
3048 'castep_keywords.json in the current directory %s '
3049 'For system wide usage, it can be copied into an ase installation '
3050 'at ASE/calculators. '
3051 'This castep_keywords.json usually only needs to be generated once '
3052 'for a CASTEP binary/CASTEP version.' % os.getcwd())
3054 import optparse
3055 parser = optparse.OptionParser()
3056 parser.add_option(
3057 '-f', '--force-write', dest='force_write',
3058 help='Force overwriting existing castep_keywords.json', default=False,
3059 action='store_true')
3060 (options, args) = parser.parse_args()
3062 if args:
3063 opt_castep_command = ''.join(args)
3064 else:
3065 opt_castep_command = ''
3066 generated = create_castep_keywords(get_castep_command(opt_castep_command),
3067 force_write=options.force_write)
3069 if generated:
3070 try:
3071 with open('castep_keywords.json') as fd:
3072 json.load(fd)
3073 except Exception as e:
3074 warnings.warn(
3075 f'{e} Ooops, something went wrong with the CASTEP keywords')
3076 else:
3077 warnings.warn('Import works. Looking good!')