Coverage for /builds/ase/ase/ase/calculators/castep.py : 41.20%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""This module defines an interface to CASTEP for
2 use by the ASE (Webpage: http://wiki.fysik.dtu.dk/ase)
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 numpy as np
18import os
19import re
20import glob
21import shutil
22import sys
23import json
24import time
25import tempfile
26import warnings
27import subprocess
28from copy import deepcopy
29from collections import namedtuple
30from itertools import product
31from pathlib import Path
32from typing import List, Set
34import ase
35import ase.units as units
36from ase.calculators.general import Calculator
37from ase.calculators.calculator import compare_atoms
38from ase.calculators.calculator import PropertyNotImplementedError
39from ase.calculators.calculator import kpts2sizeandoffsets
40from ase.dft.kpoints import BandPath
41from ase.parallel import paropen
42from ase.io.castep import read_param
43from ase.io.castep import read_bands
44from ase.constraints import FixCartesian
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 += ' {0} {1:.3f} {2:.3f} {3:.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('Cannot find band file "{}".'.format(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()
740 pass
742 # Case 2: list of explicit k-points with weights
743 # e.g. [[ 0, 0, 0, 0.125],
744 # [ 0, -0.5, 0, 0.375],
745 # [-0.5, 0, -0.5, 0.375],
746 # [-0.5, -0.5, -0.5, 0.125]]
748 elif (isinstance(kpts, (tuple, list))
749 and isinstance(kpts[0], (tuple, list))):
751 if not all(map((lambda row: len(row) == 4), kpts)):
752 raise ValueError(
753 'In explicit kpt list each row should have 4 elements')
755 clear_mp_keywords()
756 self.cell.kpoint_list = [' '.join(map(str, row)) for row in kpts]
758 # Case 3: list of explicit kpts formatted as list of str
759 # i.e. the internal format of calc.kpoint_list split on \n
760 # e.g. ['0 0 0 0.125', '0 -0.5 0 0.375', '-0.5 0 -0.5 0.375']
761 elif isinstance(kpts, (tuple, list)) and isinstance(kpts[0], str):
763 if not all(map((lambda row: len(row.split()) == 4), kpts)):
764 raise ValueError(
765 'In explicit kpt list each row should have 4 elements')
767 clear_mp_keywords()
768 self.cell.kpoint_list = kpts
770 # Case 4: list or tuple of MP samples e.g. [3, 3, 2]
771 elif isinstance(kpts, (tuple, list)) and isinstance(kpts[0], int):
772 if len(kpts) != 3:
773 raise ValueError('Monkhorst-pack grid should have 3 values')
774 clear_mp_keywords()
775 self.cell.kpoint_mp_grid = '%d %d %d' % tuple(kpts)
777 # Case 5: str representation of Case 3 e.g. '3 3 2'
778 elif isinstance(kpts, str):
779 self.set_kpts([int(x) for x in kpts.split()])
781 # Case 6: dict of options e.g. {'size': (3, 3, 2), 'gamma': True}
782 # 'spacing' is allowed but transformed to 'density' to get mesh/offset
783 elif isinstance(kpts, dict):
784 kpts = kpts.copy()
786 if (kpts.get('spacing') is not None
787 and kpts.get('density') is not None):
788 raise ValueError(
789 'Cannot set kpts spacing and density simultaneously.')
790 else:
791 if kpts.get('spacing') is not None:
792 kpts = kpts.copy()
793 spacing = kpts.pop('spacing')
794 kpts['density'] = 1 / (np.pi * spacing)
796 clear_mp_keywords()
797 size, offsets = kpts2sizeandoffsets(atoms=self.atoms, **kpts)
798 self.cell.kpoint_mp_grid = '%d %d %d' % tuple(size)
799 self.cell.kpoint_mp_offset = '%f %f %f' % tuple(offsets)
801 # Case 7: some other iterator. Try treating as a list:
802 elif hasattr(kpts, '__iter__'):
803 self.set_kpts(list(kpts))
805 # Otherwise, give up
806 else:
807 raise TypeError('Cannot interpret kpts of this type')
809 def todict(self, skip_default=True):
810 """Create dict with settings of .param and .cell"""
811 dct = {}
812 dct['param'] = self.param.get_attr_dict()
813 dct['cell'] = self.cell.get_attr_dict()
815 return dct
817 def check_state(self, atoms, tol=1e-15):
818 """Check for system changes since last calculation."""
819 return compare_atoms(self._old_atoms, atoms)
821 def _castep_find_last_record(self, castep_file):
822 """Checks wether a given castep file has a regular
823 ending message following the last banner message. If this
824 is the case, the line number of the last banner is message
825 is return, otherwise False.
827 returns (record_start, record_end, end_found, last_record_complete)
828 """
829 if isinstance(castep_file, str):
830 castep_file = paropen(castep_file, 'r')
831 file_opened = True
832 else:
833 file_opened = False
834 record_starts = []
835 while True:
836 line = castep_file.readline()
837 if (('Welcome' in line or 'Materials Studio' in line)
838 and 'CASTEP' in line):
839 record_starts = [castep_file.tell()] + record_starts
840 if not line:
841 break
843 if record_starts == []:
844 warnings.warn('Could not find CASTEP label in result file: %s.'
845 ' Are you sure this is a .castep file?' % castep_file)
846 return
848 # search for regular end of file
849 end_found = False
850 # start to search from record beginning from the back
851 # and see if
852 record_end = -1
853 for record_nr, record_start in enumerate(record_starts):
854 castep_file.seek(record_start)
855 while True:
856 line = castep_file.readline()
857 if not line:
858 break
859 if 'warn' in line.lower():
860 self._warnings.append(line)
862 if 'Finalisation time =' in line:
863 end_found = True
864 record_end = castep_file.tell()
865 break
867 if end_found:
868 break
870 if file_opened:
871 castep_file.close()
873 if end_found:
874 # record_nr == 0 corresponds to the last record here
875 if record_nr == 0:
876 return (record_start, record_end, True, True)
877 else:
878 return (record_start, record_end, True, False)
879 else:
880 return (0, record_end, False, False)
882 def read(self, castep_file=None):
883 """Read a castep file into the current instance."""
885 _close = True
887 if castep_file is None:
888 if self._castep_file:
889 castep_file = self._castep_file
890 out = paropen(castep_file, 'r')
891 else:
892 warnings.warn('No CASTEP file specified')
893 return
894 if not os.path.exists(castep_file):
895 warnings.warn('No CASTEP file found')
897 elif isinstance(castep_file, str):
898 out = paropen(castep_file, 'r')
900 else:
901 # in this case we assume that we have a fileobj already, but check
902 # for attributes in order to avoid extended EAFP blocks.
903 out = castep_file
905 # look before you leap...
906 attributes = ['name',
907 'seek',
908 'close',
909 'readline',
910 'tell']
912 for attr in attributes:
913 if not hasattr(out, attr):
914 raise TypeError(
915 '"castep_file" is neither str nor valid fileobj')
917 castep_file = out.name
918 _close = False
920 if self._seed is None:
921 self._seed = os.path.splitext(os.path.basename(castep_file))[0]
923 err_file = '%s.0001.err' % self._seed
924 if os.path.exists(err_file):
925 err_file = paropen(err_file)
926 self._error = err_file.read()
927 err_file.close()
928 # we return right-away because it might
929 # just be here from a previous run
930 # look for last result, if several CASTEP
931 # run are appended
933 record_start, record_end, end_found, _\
934 = self._castep_find_last_record(out)
935 if not end_found:
936 warnings.warn('No regular end found in %s file. %s' %
937 (castep_file, self._error))
938 if _close:
939 out.close()
940 return
941 # we return here, because the file has no a regular end
943 # now iterate over last CASTEP output in file to extract information
944 # could be generalized as well to extract trajectory from file
945 # holding several outputs
946 n_cell_const = 0
947 forces = []
949 # HOTFIX:
950 # we have to initialize the _stress variable as a zero array
951 # otherwise the calculator crashes upon pickling trajectories
952 # Alternative would be to raise a NotImplementedError() which
953 # is also kind of not true, since we can extract stresses if
954 # the user configures CASTEP to print them in the outfile
955 # stress = []
956 stress = np.zeros([3, 3])
957 hirsh_volrat = []
959 # Two flags to check whether spin-polarized or not, and whether
960 # Hirshfeld volumes are calculated
961 spin_polarized = False
962 calculate_hirshfeld = False
963 mulliken_analysis = False
964 hirshfeld_analysis = False
965 kpoints = None
967 positions_frac_list = []
968 mulliken_charges = []
969 spins = []
971 out.seek(record_start)
972 while True:
973 # TODO: add a switch if we have a geometry optimization: record
974 # atoms objects for intermediate steps.
975 try:
976 line = out.readline()
977 if not line or out.tell() > record_end:
978 break
979 elif 'Hirshfeld Analysis' in line:
980 hirshfeld_charges = []
982 hirshfeld_analysis = True
983 # skip the separating line
984 line = out.readline()
985 # this is the headline
986 line = out.readline()
988 if 'Charge' in line:
989 # skip the next separator line
990 line = out.readline()
991 while True:
992 line = out.readline()
993 fields = line.split()
994 if len(fields) == 1:
995 break
996 else:
997 hirshfeld_charges.append(float(fields[-1]))
998 elif 'stress calculation' in line:
999 if line.split()[-1].strip() == 'on':
1000 self.param.calculate_stress = True
1001 elif 'basis set accuracy' in line:
1002 self.param.basis_precision = line.split()[-1]
1003 elif 'plane wave basis set cut-off' in line:
1004 # NB this is set as a private "result" attribute to avoid
1005 # conflict with input option basis_precision
1006 cutoff = float(line.split()[-2])
1007 self._cut_off_energy = cutoff
1008 if self.param.basis_precision.value is None:
1009 self.param.cut_off_energy = cutoff
1010 elif 'total energy / atom convergence tol.' in line:
1011 elec_energy_tol = float(line.split()[-2])
1012 self.param.elec_energy_tol = elec_energy_tol
1013 elif 'convergence tolerance window' in line:
1014 elec_convergence_win = int(line.split()[-2])
1015 self.param.elec_convergence_win = elec_convergence_win
1016 elif re.match(r'\sfinite basis set correction\s*:', line):
1017 finite_basis_corr = line.split()[-1]
1018 fbc_possibilities = {'none': 0,
1019 'manual': 1, 'automatic': 2}
1020 fbc = fbc_possibilities[finite_basis_corr]
1021 self.param.finite_basis_corr = fbc
1022 elif 'Treating system as non-metallic' in line:
1023 self.param.fix_occupancy = True
1024 elif 'max. number of SCF cycles:' in line:
1025 max_no_scf = float(line.split()[-1])
1026 self.param.max_scf_cycles = max_no_scf
1027 elif 'density-mixing scheme' in line:
1028 mixing_scheme = line.split()[-1]
1029 self.param.mixing_scheme = mixing_scheme
1030 elif 'dump wavefunctions every' in line:
1031 no_dump_cycles = float(line.split()[-3])
1032 self.param.num_dump_cycles = no_dump_cycles
1033 elif 'optimization strategy' in line:
1034 lspl = line.split(":")
1035 if lspl[0].strip() != 'optimization strategy':
1036 # This can happen in iprint: 3 calculations
1037 continue
1038 if 'memory' in line:
1039 self.param.opt_strategy = 'Memory'
1040 if 'speed' in line:
1041 self.param.opt_strategy = 'Speed'
1042 elif 'calculation limited to maximum' in line:
1043 calc_limit = float(line.split()[-2])
1044 self.param.run_time = calc_limit
1045 elif 'type of calculation' in line:
1046 lspl = line.split(":")
1047 if lspl[0].strip() != 'type of calculation':
1048 # This can happen in iprint: 3 calculations
1049 continue
1050 calc_type = lspl[-1]
1051 calc_type = re.sub(r'\s+', ' ', calc_type)
1052 calc_type = calc_type.strip()
1053 if calc_type != 'single point energy':
1054 calc_type_possibilities = {
1055 'geometry optimization': 'GeometryOptimization',
1056 'band structure': 'BandStructure',
1057 'molecular dynamics': 'MolecularDynamics',
1058 'optical properties': 'Optics',
1059 'phonon calculation': 'Phonon',
1060 'E-field calculation': 'Efield',
1061 'Phonon followed by E-field': 'Phonon+Efield',
1062 'transition state search': 'TransitionStateSearch',
1063 'Magnetic Resonance': 'MagRes',
1064 'Core level spectra': 'Elnes',
1065 'Electronic Spectroscopy': 'ElectronicSpectroscopy'
1066 }
1067 ctype = calc_type_possibilities[calc_type]
1068 self.param.task = ctype
1069 elif 'using functional' in line:
1070 used_functional = line.split(":")[-1]
1071 used_functional = re.sub(r'\s+', ' ', used_functional)
1072 used_functional = used_functional.strip()
1073 if used_functional != 'Local Density Approximation':
1074 used_functional_possibilities = {
1075 'Perdew Wang (1991)': 'PW91',
1076 'Perdew Burke Ernzerhof': 'PBE',
1077 'revised Perdew Burke Ernzerhof': 'RPBE',
1078 'PBE with Wu-Cohen exchange': 'WC',
1079 'PBE for solids (2008)': 'PBESOL',
1080 'Hartree-Fock': 'HF',
1081 'Hartree-Fock +': 'HF-LDA',
1082 'Screened Hartree-Fock': 'sX',
1083 'Screened Hartree-Fock + ': 'sX-LDA',
1084 'hybrid PBE0': 'PBE0',
1085 'hybrid B3LYP': 'B3LYP',
1086 'hybrid HSE03': 'HSE03',
1087 'hybrid HSE06': 'HSE06'
1088 }
1089 used_func = used_functional_possibilities[
1090 used_functional]
1091 self.param.xc_functional = used_func
1092 elif 'output verbosity' in line:
1093 iprint = int(line.split()[-1][1])
1094 if int(iprint) != 1:
1095 self.param.iprint = iprint
1096 elif 'treating system as spin-polarized' in line:
1097 spin_polarized = True
1098 self.param.spin_polarized = spin_polarized
1099 elif 'treating system as non-spin-polarized' in line:
1100 spin_polarized = False
1101 elif 'Number of kpoints used' in line:
1102 kpoints = int(line.split('=')[-1].strip())
1103 elif 'Unit Cell' in line:
1104 lattice_real = []
1105 lattice_reci = []
1106 while True:
1107 line = out.readline()
1108 fields = line.split()
1109 if len(fields) == 6:
1110 break
1111 for i in range(3):
1112 lattice_real.append([float(f) for f in fields[0:3]])
1113 lattice_reci.append([float(f) for f in fields[3:7]])
1114 line = out.readline()
1115 fields = line.split()
1116 elif 'Cell Contents' in line:
1117 while True:
1118 line = out.readline()
1119 if 'Total number of ions in cell' in line:
1120 n_atoms = int(line.split()[7])
1121 if 'Total number of species in cell' in line:
1122 int(line.split()[7])
1123 fields = line.split()
1124 if len(fields) == 0:
1125 break
1126 elif 'Fractional coordinates of atoms' in line:
1127 species = []
1128 custom_species = None # A CASTEP special thing
1129 positions_frac = []
1130 # positions_cart = []
1131 while True:
1132 line = out.readline()
1133 fields = line.split()
1134 if len(fields) == 7:
1135 break
1136 for n in range(n_atoms):
1137 spec_custom = fields[1].split(':', 1)
1138 elem = spec_custom[0]
1139 if len(spec_custom) > 1 and custom_species is None:
1140 # Add it to the custom info!
1141 custom_species = list(species)
1142 species.append(elem)
1143 if custom_species is not None:
1144 custom_species.append(fields[1])
1145 positions_frac.append([float(s) for s in fields[3:6]])
1146 line = out.readline()
1147 fields = line.split()
1148 positions_frac_list.append(positions_frac)
1149 elif 'Files used for pseudopotentials' in line:
1150 while True:
1151 line = out.readline()
1152 if 'Pseudopotential generated on-the-fly' in line:
1153 continue
1154 fields = line.split()
1155 if (len(fields) >= 2):
1156 elem, pp_file = fields
1157 self.cell.species_pot = (elem, pp_file)
1158 else:
1159 break
1160 elif 'k-Points For BZ Sampling' in line:
1161 # TODO: generalize for non-Monkhorst Pack case
1162 # (i.e. kpoint lists) -
1163 # kpoints_offset cannot be read this way and
1164 # is hence always set to None
1165 while True:
1166 line = out.readline()
1167 if not line.strip():
1168 break
1169 if 'MP grid size for SCF calculation' in line:
1170 # kpoints = ' '.join(line.split()[-3:])
1171 # self.kpoints_mp_grid = kpoints
1172 # self.kpoints_mp_offset = '0. 0. 0.'
1173 # not set here anymore because otherwise
1174 # two calculator objects go out of sync
1175 # after each calculation triggering unnecessary
1176 # recalculation
1177 break
1178 elif 'Number of cell constraints' in line:
1179 n_cell_const = int(line.split()[4])
1180 elif 'Final energy' in line:
1181 self._energy_total = float(line.split()[-2])
1182 elif 'Final free energy' in line:
1183 self._energy_free = float(line.split()[-2])
1184 elif 'NB est. 0K energy' in line:
1185 self._energy_0K = float(line.split()[-2])
1186 # check if we had a finite basis set correction
1187 elif 'Total energy corrected for finite basis set' in line:
1188 self._energy_total_corr = float(line.split()[-2])
1190 # Add support for dispersion correction
1191 # filtering due to SEDC is done in get_potential_energy
1192 elif 'Dispersion corrected final energy' in line:
1193 self._dispcorr_energy_total = float(line.split()[-2])
1194 elif 'Dispersion corrected final free energy' in line:
1195 self._dispcorr_energy_free = float(line.split()[-2])
1196 elif 'dispersion corrected est. 0K energy' in line:
1197 self._dispcorr_energy_0K = float(line.split()[-2])
1199 # remember to remove constraint labels in force components
1200 # (lacking a space behind the actual floating point number in
1201 # the CASTEP output)
1202 elif '******************** Forces *********************'\
1203 in line or\
1204 '************** Symmetrised Forces ***************'\
1205 in line or\
1206 '************** Constrained Symmetrised Forces ****'\
1207 '**********'\
1208 in line or\
1209 '******************** Constrained Forces **********'\
1210 '**********'\
1211 in line or\
1212 '******************* Unconstrained Forces *********'\
1213 '**********'\
1214 in line:
1215 fix = []
1216 fix_cart = []
1217 forces = []
1218 while True:
1219 line = out.readline()
1220 fields = line.split()
1221 if len(fields) == 7:
1222 break
1223 for n in range(n_atoms):
1224 consd = np.array([0, 0, 0])
1225 fxyz = [0, 0, 0]
1226 for (i, force_component) in enumerate(fields[-4:-1]):
1227 if force_component.count("(cons'd)") > 0:
1228 consd[i] = 1
1229 fxyz[i] = float(force_component.replace(
1230 "(cons'd)", ''))
1231 if consd.all():
1232 fix.append(n)
1233 elif consd.any():
1234 fix_cart.append(FixCartesian(n, consd))
1235 forces.append(fxyz)
1236 line = out.readline()
1237 fields = line.split()
1239 # add support for Hirshfeld analysis
1240 elif 'Hirshfeld / free atomic volume :' in line:
1241 # if we are here, then params must be able to cope with
1242 # Hirshfeld flag (if castep_keywords.py matches employed
1243 # castep version)
1244 calculate_hirshfeld = True
1245 hirsh_volrat = []
1246 while True:
1247 line = out.readline()
1248 fields = line.split()
1249 if len(fields) == 1:
1250 break
1251 for n in range(n_atoms):
1252 hirsh_atom = float(fields[0])
1253 hirsh_volrat.append(hirsh_atom)
1254 while True:
1255 line = out.readline()
1256 if 'Hirshfeld / free atomic volume :' in line or\
1257 'Hirshfeld Analysis' in line:
1258 break
1259 line = out.readline()
1260 fields = line.split()
1262 elif '***************** Stress Tensor *****************'\
1263 in line or\
1264 '*********** Symmetrised Stress Tensor ***********'\
1265 in line:
1266 stress = []
1267 while True:
1268 line = out.readline()
1269 fields = line.split()
1270 if len(fields) == 6:
1271 break
1272 for n in range(3):
1273 stress.append([float(s) for s in fields[2:5]])
1274 line = out.readline()
1275 fields = line.split()
1276 line = out.readline()
1277 if "Pressure:" in line:
1278 self._pressure = float(line.split()[-2]) * units.GPa
1279 elif ('BFGS: starting iteration' in line
1280 or 'BFGS: improving iteration' in line):
1281 if n_cell_const < 6:
1282 lattice_real = []
1283 lattice_reci = []
1284 # backup previous configuration first:
1285 # for highly symmetric systems (where essentially only the
1286 # stress is optimized, but the atomic positions) positions
1287 # are only printed once.
1288 if species:
1289 prev_species = deepcopy(species)
1290 if positions_frac:
1291 prev_positions_frac = deepcopy(positions_frac)
1292 species = []
1293 positions_frac = []
1294 forces = []
1296 # HOTFIX:
1297 # Same reason for the stress initialization as before
1298 # stress = []
1299 stress = np.zeros([3, 3])
1301 # extract info from the Mulliken analysis
1302 elif 'Atomic Populations' in line:
1303 # sometimes this appears twice in a castep file
1305 mulliken_analysis = True
1306 # skip the separating line
1307 line = out.readline()
1308 # this is the headline
1309 line = out.readline()
1311 if 'Charge' in line:
1312 # skip the next separator line
1313 line = out.readline()
1314 while True:
1315 line = out.readline()
1316 fields = line.split()
1317 if len(fields) == 1:
1318 break
1320 # the check for len==7 is due to CASTEP 18
1321 # outformat changes
1322 if spin_polarized:
1323 if len(fields) != 7:
1324 spins.append(float(fields[-1]))
1325 mulliken_charges.append(float(fields[-2]))
1326 else:
1327 mulliken_charges.append(float(fields[-1]))
1329 # There is actually no good reason to get out of the loop
1330 # already at this point... or do I miss something?
1331 # elif 'BFGS: Final Configuration:' in line:
1332 # break
1333 elif 'warn' in line.lower():
1334 self._warnings.append(line)
1336 # fetch some last info
1337 elif 'Total time' in line:
1338 pattern = r'.*=\s*([\d\.]+) s'
1339 self._total_time = float(re.search(pattern, line).group(1))
1341 elif 'Peak Memory Use' in line:
1342 pattern = r'.*=\s*([\d]+) kB'
1343 self._peak_memory = int(re.search(pattern, line).group(1))
1345 except Exception as exception:
1346 sys.stderr.write(line + '|-> line triggered exception: '
1347 + str(exception))
1348 raise
1350 if _close:
1351 out.close()
1353 # in highly summetric crystals, positions and symmetry are only printed
1354 # upon init, hence we here restore these original values
1355 if not positions_frac:
1356 positions_frac = prev_positions_frac
1357 if not species:
1358 species = prev_species
1360 if not spin_polarized:
1361 # set to zero spin if non-spin polarized calculation
1362 spins = np.zeros(len(positions_frac))
1363 elif len(spins) != len(positions_frac):
1364 warnings.warn('Spins could not be read for the atoms despite'
1365 ' spin-polarized calculation; spins will be ignored')
1366 spins = np.zeros(len(positions_frac))
1368 positions_frac_atoms = np.array(positions_frac)
1369 forces_atoms = np.array(forces)
1370 spins_atoms = np.array(spins)
1372 if mulliken_analysis:
1373 if len(mulliken_charges) != len(positions_frac):
1374 warnings.warn(
1375 'Mulliken charges could not be read for the atoms;'
1376 ' charges will be ignored')
1377 mulliken_charges = np.zeros(len(positions_frac))
1378 mulliken_charges_atoms = np.array(mulliken_charges)
1379 else:
1380 mulliken_charges_atoms = np.zeros(len(positions_frac))
1382 if hirshfeld_analysis:
1383 hirshfeld_charges_atoms = np.array(hirshfeld_charges)
1384 else:
1385 hirshfeld_charges_atoms = None
1387 if calculate_hirshfeld:
1388 hirsh_atoms = np.array(hirsh_volrat)
1389 else:
1390 hirsh_atoms = np.zeros_like(spins)
1392 if self.atoms and not self._set_atoms:
1393 # compensate for internal reordering of atoms by CASTEP
1394 # using the fact that the order is kept within each species
1396 # positions_frac_ase = self.atoms.get_scaled_positions(wrap=False)
1397 atoms_assigned = [False] * len(self.atoms)
1399 # positions_frac_castep_init = np.array(positions_frac_list[0])
1400 positions_frac_castep = np.array(positions_frac_list[-1])
1402 # species_castep = list(species)
1403 forces_castep = np.array(forces)
1404 hirsh_castep = np.array(hirsh_volrat)
1405 spins_castep = np.array(spins)
1406 mulliken_charges_castep = np.array(mulliken_charges_atoms)
1408 # go through the atoms position list and replace
1409 # with the corresponding one from the
1410 # castep file corresponding atomic number
1411 for iase in range(n_atoms):
1412 for icastep in range(n_atoms):
1413 if (species[icastep] == self.atoms[iase].symbol
1414 and not atoms_assigned[icastep]):
1415 positions_frac_atoms[iase] = \
1416 positions_frac_castep[icastep]
1417 forces_atoms[iase] = np.array(forces_castep[icastep])
1418 if iprint > 1 and calculate_hirshfeld:
1419 hirsh_atoms[iase] = np.array(hirsh_castep[icastep])
1420 if spin_polarized:
1421 # reordering not necessary in case all spins == 0
1422 spins_atoms[iase] = np.array(spins_castep[icastep])
1423 mulliken_charges_atoms[iase] = np.array(
1424 mulliken_charges_castep[icastep])
1425 atoms_assigned[icastep] = True
1426 break
1428 if not all(atoms_assigned):
1429 not_assigned = [i for (i, assigned)
1430 in zip(range(len(atoms_assigned)),
1431 atoms_assigned) if not assigned]
1432 warnings.warn(
1433 '%s atoms not assigned. '
1434 ' DEBUGINFO: The following atoms where not assigned: %s' %
1435 (atoms_assigned.count(False), not_assigned))
1436 else:
1437 self.atoms.set_scaled_positions(positions_frac_atoms)
1439 else:
1440 # If no atoms, object has been previously defined
1441 # we define it here and set the Castep() instance as calculator.
1442 # This covers the case that we simply want to open a .castep file.
1444 # The next time around we will have an atoms object, since
1445 # set_calculator also set atoms in the calculator.
1446 if self.atoms:
1447 constraints = self.atoms.constraints
1448 else:
1449 constraints = []
1450 atoms = ase.atoms.Atoms(species,
1451 cell=lattice_real,
1452 constraint=constraints,
1453 pbc=True,
1454 scaled_positions=positions_frac,
1455 )
1456 if custom_species is not None:
1457 atoms.new_array('castep_custom_species',
1458 np.array(custom_species))
1460 if self.param.spin_polarized:
1461 # only set magnetic moments if this was a spin polarized
1462 # calculation
1463 # this one fails as is
1464 atoms.set_initial_magnetic_moments(magmoms=spins_atoms)
1466 if mulliken_analysis:
1467 atoms.set_initial_charges(charges=mulliken_charges_atoms)
1468 atoms.calc = self
1470 self._kpoints = kpoints
1471 self._forces = forces_atoms
1472 # stress in .castep file is given in GPa:
1473 self._stress = np.array(stress) * units.GPa
1474 self._hirsh_volrat = hirsh_atoms
1475 self._spins = spins_atoms
1476 self._mulliken_charges = mulliken_charges_atoms
1477 self._hirshfeld_charges = hirshfeld_charges_atoms
1479 if self._warnings:
1480 warnings.warn('WARNING: %s contains warnings' % castep_file)
1481 for warning in self._warnings:
1482 warnings.warn(warning)
1483 # reset
1484 self._warnings = []
1486 # Read in eigenvalues from bands file
1487 bands_file = castep_file[:-7] + '.bands'
1488 if (self.param.task.value is not None
1489 and self.param.task.value.lower() == 'bandstructure'):
1490 self._band_structure = self.band_structure(bandfile=bands_file)
1491 else:
1492 try:
1493 (self._ibz_kpts,
1494 self._ibz_weights,
1495 self._eigenvalues,
1496 self._efermi) = read_bands(filename=bands_file)
1497 except FileNotFoundError:
1498 warnings.warn('Could not load .bands file, eigenvalues and '
1499 'Fermi energy are unknown')
1501 def get_hirsh_volrat(self):
1502 """
1503 Return the Hirshfeld volumes.
1504 """
1505 return self._hirsh_volrat
1507 def get_spins(self):
1508 """
1509 Return the spins from a plane-wave Mulliken analysis.
1510 """
1511 return self._spins
1513 def get_mulliken_charges(self):
1514 """
1515 Return the charges from a plane-wave Mulliken analysis.
1516 """
1517 return self._mulliken_charges
1519 def get_hirshfeld_charges(self):
1520 """
1521 Return the charges from a Hirshfeld analysis.
1522 """
1523 return self._hirshfeld_charges
1525 def get_total_time(self):
1526 """
1527 Return the total runtime
1528 """
1529 return self._total_time
1531 def get_peak_memory(self):
1532 """
1533 Return the peak memory usage
1534 """
1535 return self._peak_memory
1537 def set_label(self, label):
1538 """The label is part of each seed, which in turn is a prefix
1539 in each CASTEP related file.
1540 """
1541 # we may think about changing this in future to set `self._directory`
1542 # and `self._label`, as one would expect
1543 self._label = label
1545 def set_pspot(self, pspot, elems=None,
1546 notelems=None,
1547 clear=True,
1548 suffix='usp'):
1549 """Quickly set all pseudo-potentials: Usually CASTEP psp are named
1550 like <Elem>_<pspot>.<suffix> so this function function only expects
1551 the <LibraryName>. It then clears any previous pseudopotential
1552 settings apply the one with <LibraryName> for each element in the
1553 atoms object. The optional elems and notelems arguments can be used
1554 to exclusively assign to some species, or to exclude with notelemens.
1556 Parameters ::
1558 - elems (None) : set only these elements
1559 - notelems (None): do not set the elements
1560 - clear (True): clear previous settings
1561 - suffix (usp): PP file suffix
1562 """
1563 if self._find_pspots:
1564 if self._pedantic:
1565 warnings.warn('Warning: <_find_pspots> = True. '
1566 'Do you really want to use `set_pspots()`? '
1567 'This does not check whether the PP files exist. '
1568 'You may rather want to use `find_pspots()` with '
1569 'the same <pspot>.')
1571 if clear and not elems and not notelems:
1572 self.cell.species_pot.clear()
1573 for elem in set(self.atoms.get_chemical_symbols()):
1574 if elems is not None and elem not in elems:
1575 continue
1576 if notelems is not None and elem in notelems:
1577 continue
1578 self.cell.species_pot = (elem, '%s_%s.%s' % (elem, pspot, suffix))
1580 def find_pspots(self, pspot='.+', elems=None,
1581 notelems=None, clear=True, suffix='(usp|UPF|recpot)'):
1582 r"""Quickly find and set all pseudo-potentials by searching in
1583 castep_pp_path:
1585 This one is more flexible than set_pspots, and also checks if the files
1586 are actually available from the castep_pp_path.
1588 Essentially, the function parses the filenames in <castep_pp_path> and
1589 does a regex matching. The respective pattern is:
1591 r"^(<elem>|<elem.upper()>|elem.lower()>(_|-)<pspot>\.<suffix>$"
1593 In most cases, it will be sufficient to not specify anything, if you
1594 use standard CASTEP USPPs with only one file per element in the
1595 <castep_pp_path>.
1597 The function raises a `RuntimeError` if there is some ambiguity
1598 (multiple files per element).
1600 Parameters ::
1602 - pspots ('.+') : as defined above, will be a wildcard if not
1603 specified.
1604 - elems (None) : set only these elements
1605 - notelems (None): do not set the elements
1606 - clear (True): clear previous settings
1607 - suffix (usp|UPF|recpot): PP file suffix
1608 """
1609 if clear and not elems and not notelems:
1610 self.cell.species_pot.clear()
1612 if not os.path.isdir(self._castep_pp_path):
1613 if self._pedantic:
1614 warnings.warn(
1615 'Cannot search directory: {} Folder does not exist'
1616 .format(self._castep_pp_path))
1617 return
1619 # translate the bash wildcard syntax to regex
1620 if pspot == '*':
1621 pspot = '.*'
1622 if suffix == '*':
1623 suffix = '.*'
1624 if pspot == '*':
1625 pspot = '.*'
1627 # GBRV USPPs have a strnage naming schme
1628 pattern = r'^({elem}|{elem_upper}|{elem_lower})(_|-){pspot}\.{suffix}$'
1630 for elem in set(self.atoms.get_chemical_symbols()):
1631 if elems is not None and elem not in elems:
1632 continue
1633 if notelems is not None and elem in notelems:
1634 continue
1635 p = pattern.format(elem=elem,
1636 elem_upper=elem.upper(),
1637 elem_lower=elem.lower(),
1638 pspot=pspot,
1639 suffix=suffix)
1640 pps = []
1641 for f in os.listdir(self._castep_pp_path):
1642 if re.match(p, f):
1643 pps.append(f)
1644 if not pps:
1645 if self._pedantic:
1646 warnings.warn('Pseudopotential for species {} not found!'
1647 .format(elem))
1648 elif not len(pps) == 1:
1649 raise RuntimeError(
1650 'Pseudopotential for species ''{} not unique!\n'
1651 .format(elem)
1652 + 'Found the following files in {}\n'
1653 .format(self._castep_pp_path)
1654 + '\n'.join([' {}'.format(pp) for pp in pps]) +
1655 '\nConsider a stricter search pattern in `find_pspots()`.')
1656 else:
1657 self.cell.species_pot = (elem, pps[0])
1659 @property
1660 def name(self):
1661 """Return the name of the calculator (string). """
1662 return self.__name__
1664 def get_property(self, name, atoms=None, allow_calculation=True):
1665 # High-level getter for compliance with the database module...
1666 # in principle this would not be necessary any longer if we properly
1667 # based this class on `Calculator`
1668 if name == 'forces':
1669 return self.get_forces(atoms)
1670 elif name == 'energy':
1671 return self.get_potential_energy(atoms)
1672 elif name == 'stress':
1673 return self.get_stress(atoms)
1674 elif name == 'charges':
1675 return self.get_charges(atoms)
1676 else:
1677 raise PropertyNotImplementedError
1679 @_self_getter
1680 def get_forces(self, atoms):
1681 """Run CASTEP calculation if needed and return forces."""
1682 self.update(atoms)
1683 return np.array(self._forces)
1685 @_self_getter
1686 def get_total_energy(self, atoms):
1687 """Run CASTEP calculation if needed and return total energy."""
1688 self.update(atoms)
1689 return self._energy_total
1691 @_self_getter
1692 def get_total_energy_corrected(self, atoms):
1693 """Run CASTEP calculation if needed and return total energy."""
1694 self.update(atoms)
1695 return self._energy_total_corr
1697 @_self_getter
1698 def get_free_energy(self, atoms):
1699 """Run CASTEP calculation if needed and return free energy.
1700 Only defined with smearing."""
1701 self.update(atoms)
1702 return self._energy_free
1704 @_self_getter
1705 def get_0K_energy(self, atoms):
1706 """Run CASTEP calculation if needed and return 0K energy.
1707 Only defined with smearing."""
1708 self.update(atoms)
1709 return self._energy_0K
1711 @_self_getter
1712 def get_potential_energy(self, atoms, force_consistent=False):
1713 # here for compatibility with ase/calculators/general.py
1714 # but accessing only _name variables
1715 """Return the total potential energy."""
1716 self.update(atoms)
1717 if force_consistent:
1718 # Assumption: If no dispersion correction is applied, then the
1719 # respective value will default to None as initialized.
1720 if self._dispcorr_energy_free is not None:
1721 return self._dispcorr_energy_free
1722 else:
1723 return self._energy_free
1724 else:
1725 if self._energy_0K is not None:
1726 if self._dispcorr_energy_0K is not None:
1727 return self._dispcorr_energy_0K
1728 else:
1729 return self._energy_0K
1730 else:
1731 if self._dispcorr_energy_total is not None:
1732 return self._dispcorr_energy_total
1733 else:
1734 if self._energy_total_corr is not None:
1735 return self._energy_total_corr
1736 else:
1737 return self._energy_total
1739 @_self_getter
1740 def get_stress(self, atoms):
1741 """Return the stress."""
1742 self.update(atoms)
1743 # modification: we return the Voigt form directly to get rid of the
1744 # annoying user warnings
1745 stress = np.array(
1746 [self._stress[0, 0], self._stress[1, 1], self._stress[2, 2],
1747 self._stress[1, 2], self._stress[0, 2], self._stress[0, 1]])
1748 # return self._stress
1749 return stress
1751 @_self_getter
1752 def get_pressure(self, atoms):
1753 """Return the pressure."""
1754 self.update(atoms)
1755 return self._pressure
1757 @_self_getter
1758 def get_unit_cell(self, atoms):
1759 """Return the unit cell."""
1760 self.update(atoms)
1761 return self._unit_cell
1763 @_self_getter
1764 def get_kpoints(self, atoms):
1765 """Return the kpoints."""
1766 self.update(atoms)
1767 return self._kpoints
1769 @_self_getter
1770 def get_number_cell_constraints(self, atoms):
1771 """Return the number of cell constraints."""
1772 self.update(atoms)
1773 return self._number_of_cell_constraints
1775 @_self_getter
1776 def get_charges(self, atoms):
1777 """Run CASTEP calculation if needed and return Mulliken charges."""
1778 self.update(atoms)
1779 return np.array(self._mulliken_charges)
1781 @_self_getter
1782 def get_magnetic_moments(self, atoms):
1783 """Run CASTEP calculation if needed and return Mulliken charges."""
1784 self.update(atoms)
1785 return np.array(self._spins)
1787 def set_atoms(self, atoms):
1788 """Sets the atoms for the calculator and vice versa."""
1789 atoms.pbc = [True, True, True]
1790 self.__dict__['atoms'] = atoms.copy()
1791 self.atoms._calc = self
1793 def update(self, atoms):
1794 """Checks if atoms object or calculator changed and
1795 runs calculation if so.
1796 """
1797 if self.calculation_required(atoms):
1798 self.calculate(atoms)
1800 def calculation_required(self, atoms, _=None):
1801 """Checks wether anything changed in the atoms object or CASTEP
1802 settings since the last calculation using this instance.
1803 """
1804 # SPR: what happens with the atoms parameter here? Why don't we use it?
1805 # from all that I can tell we need to compare against atoms instead of
1806 # self.atoms
1807 # if not self.atoms == self._old_atoms:
1808 if not atoms == self._old_atoms:
1809 return True
1810 if self._old_param is None or self._old_cell is None:
1811 return True
1812 if not self.param._options == self._old_param._options:
1813 return True
1814 if not self.cell._options == self._old_cell._options:
1815 return True
1816 return False
1818 def calculate(self, atoms):
1819 """Write all necessary input file and call CASTEP."""
1820 self.prepare_input_files(atoms, force_write=self._force_write)
1821 if not self._prepare_input_only:
1822 self.run()
1823 self.read()
1825 # we need to push the old state here!
1826 # although run() pushes it, read() may change the atoms object
1827 # again.
1828 # yet, the old state is supposed to be the one AFTER read()
1829 self.push_oldstate()
1831 def push_oldstate(self):
1832 """This function pushes the current state of the (CASTEP) Atoms object
1833 onto the previous state. Or in other words after calling this function,
1834 calculation_required will return False and enquiry functions just
1835 report the current value, e.g. get_forces(), get_potential_energy().
1836 """
1837 # make a snapshot of all current input
1838 # to be able to test if recalculation
1839 # is necessary
1840 self._old_atoms = self.atoms.copy()
1841 self._old_param = deepcopy(self.param)
1842 self._old_cell = deepcopy(self.cell)
1844 def initialize(self, *args, **kwargs):
1845 """Just an alias for prepar_input_files to comply with standard
1846 function names in ASE.
1847 """
1848 self.prepare_input_files(*args, **kwargs)
1850 def prepare_input_files(self, atoms=None, force_write=None):
1851 """Only writes the input .cell and .param files and return
1852 This can be useful if one quickly needs to prepare input files
1853 for a cluster where no python or ASE is available. One can than
1854 upload the file manually and read out the results using
1855 Castep().read().
1856 """
1858 if self.param.reuse.value is None:
1859 if self._pedantic:
1860 warnings.warn(
1861 'You have not set e.g. calc.param.reuse = True. '
1862 'Reusing a previous calculation may save CPU time! '
1863 'The interface will make sure by default, .check exists. '
1864 'file before adding this statement to the .param file.')
1865 if self.param.num_dump_cycles.value is None:
1866 if self._pedantic:
1867 warnings.warn(
1868 'You have not set e.g. calc.param.num_dump_cycles = 0. '
1869 'This can save you a lot of disk space. One only needs '
1870 '*wvfn* if electronic convergence is not achieved.')
1871 from ase.io.castep import write_param
1873 if atoms is None:
1874 atoms = self.atoms
1875 else:
1876 self.atoms = atoms
1878 if force_write is None:
1879 force_write = self._force_write
1881 # if we have new instance of the calculator,
1882 # move existing results out of the way, first
1883 if (os.path.isdir(self._directory)
1884 and self._calls == 0
1885 and self._rename_existing_dir):
1886 if os.listdir(self._directory) == []:
1887 os.rmdir(self._directory)
1888 else:
1889 # rename appending creation date of the directory
1890 ctime = time.localtime(os.lstat(self._directory).st_ctime)
1891 os.rename(self._directory, '%s.bak-%s' %
1892 (self._directory,
1893 time.strftime('%Y%m%d-%H%M%S', ctime)))
1895 # create work directory
1896 if not os.path.isdir(self._directory):
1897 os.makedirs(self._directory, 0o775)
1899 # we do this every time, not only upon first call
1900 # if self._calls == 0:
1901 self._fetch_pspots()
1903 # if _try_reuse is requested and this
1904 # is not the first run, we try to find
1905 # the .check file from the previous run
1906 # this is only necessary if _track_output
1907 # is set to true
1908 if self._try_reuse and self._calls > 0:
1909 if os.path.exists(self._abs_path(self._check_file)):
1910 self.param.reuse = self._check_file
1911 elif os.path.exists(self._abs_path(self._castep_bin_file)):
1912 self.param.reuse = self._castep_bin_file
1913 self._seed = self._build_castep_seed()
1914 self._check_file = '%s.check' % self._seed
1915 self._castep_bin_file = '%s.castep_bin' % self._seed
1916 self._castep_file = self._abs_path('%s.castep' % self._seed)
1918 # write out the input file
1919 self._write_cell(self._abs_path('%s.cell' % self._seed),
1920 self.atoms, castep_cell=self.cell,
1921 force_write=force_write)
1923 if self._export_settings:
1924 interface_options = self._opt
1925 else:
1926 interface_options = None
1927 write_param(self._abs_path('%s.param' % self._seed), self.param,
1928 check_checkfile=self._check_checkfile,
1929 force_write=force_write,
1930 interface_options=interface_options,)
1932 def _build_castep_seed(self):
1933 """Abstracts to construction of the final castep <seed>
1934 with and without _tracking_output.
1935 """
1936 if self._track_output:
1937 return '%s-%06d' % (self._label, self._calls)
1938 else:
1939 return '%s' % (self._label)
1941 def _abs_path(self, path):
1942 # Create an absolute path for a file to put in the working directory
1943 return os.path.join(self._directory, path)
1945 def run(self):
1946 """Simply call castep. If the first .err file
1947 contains text, this will be printed to the screen.
1948 """
1949 # change to target directory
1950 self._calls += 1
1952 # run castep itself
1953 stdout, stderr = shell_stdouterr('%s %s' % (self._castep_command,
1954 self._seed),
1955 cwd=self._directory)
1956 if stdout:
1957 print('castep call stdout:\n%s' % stdout)
1958 if stderr:
1959 print('castep call stderr:\n%s' % stderr)
1961 # shouldn't it be called after read()???
1962 # self.push_oldstate()
1964 # check for non-empty error files
1965 err_file = self._abs_path('%s.0001.err' % self._seed)
1966 if os.path.exists(err_file):
1967 err_file = open(err_file)
1968 self._error = err_file.read()
1969 err_file.close()
1970 if self._error:
1971 raise RuntimeError(self._error)
1973 def __repr__(self):
1974 """Returns generic, fast to capture representation of
1975 CASTEP settings along with atoms object.
1976 """
1977 expr = ''
1978 expr += '-----------------Atoms--------------------\n'
1979 if self.atoms is not None:
1980 expr += str('%20s\n' % self.atoms)
1981 else:
1982 expr += 'None\n'
1984 expr += '-----------------Param keywords-----------\n'
1985 expr += str(self.param)
1986 expr += '-----------------Cell keywords------------\n'
1987 expr += str(self.cell)
1988 expr += '-----------------Internal keys------------\n'
1989 for key in self.internal_keys:
1990 expr += '%20s : %s\n' % (key, self._opt[key])
1992 return expr
1994 def __getattr__(self, attr):
1995 """___getattr___ gets overloaded to reroute the internal keys
1996 and to be able to easily store them in in the param so that
1997 they can be read in again in subsequent calls.
1998 """
1999 if attr in self.internal_keys:
2000 return self._opt[attr]
2001 if attr in ['__repr__', '__str__']:
2002 raise AttributeError
2003 elif attr not in self.__dict__:
2004 raise AttributeError('Attribute {0} not found'.format(attr))
2005 else:
2006 return self.__dict__[attr]
2008 def __setattr__(self, attr, value):
2009 """We overload the settattr method to make value assignment
2010 as pythonic as possible. Internal values all start with _.
2011 Value assigment is case insensitive!
2012 """
2014 if attr.startswith('_'):
2015 # internal variables all start with _
2016 # let's check first if they are close but not identical
2017 # to one of the switches, that the user accesses directly
2018 similars = difflib.get_close_matches(attr, self.internal_keys,
2019 cutoff=0.9)
2020 if attr not in self.internal_keys and similars:
2021 warnings.warn(
2022 'Warning: You probably tried one of: %s but typed %s' %
2023 (similars, attr))
2024 if attr in self.internal_keys:
2025 self._opt[attr] = value
2026 if attr == '_track_output':
2027 if value:
2028 self._try_reuse = True
2029 if self._pedantic:
2030 warnings.warn(
2031 'You switched _track_output on. This will '
2032 'consume a lot of disk-space. The interface '
2033 'also switched _try_reuse on, which will '
2034 'try to find the last check file. Set '
2035 '_try_reuse = False, if you need '
2036 'really separate calculations')
2037 elif '_try_reuse' in self._opt and self._try_reuse:
2038 self._try_reuse = False
2039 if self._pedantic:
2040 warnings.warn('_try_reuse is set to False, too')
2041 else:
2042 self.__dict__[attr] = value
2043 return
2044 elif attr in ['atoms', 'cell', 'param']:
2045 if value is not None:
2046 if attr == 'atoms' and not isinstance(value, ase.atoms.Atoms):
2047 raise TypeError(
2048 '%s is not an instance of ase.atoms.Atoms.' % value)
2049 elif attr == 'cell' and not isinstance(value, CastepCell):
2050 raise TypeError('%s is not an instance of CastepCell.' %
2051 value)
2052 elif attr == 'param' and not isinstance(value, CastepParam):
2053 raise TypeError('%s is not an instance of CastepParam.' %
2054 value)
2055 # These 3 are accepted right-away, no matter what
2056 self.__dict__[attr] = value
2057 return
2058 elif attr in self.atoms_obj_keys:
2059 # keywords which clearly belong to the atoms object are
2060 # rerouted to go there
2061 self.atoms.__dict__[attr] = value
2062 return
2063 elif attr in self.atoms_keys:
2064 # CASTEP keywords that should go into the atoms object
2065 # itself are blocked
2066 warnings.warn('Ignoring setings of "%s", since this has to be set '
2067 'through the atoms object' % attr)
2068 return
2070 attr = attr.lower()
2071 if attr not in (list(self.cell._options.keys())
2072 + list(self.param._options.keys())):
2073 # what is left now should be meant to be a castep keyword
2074 # so we first check if it defined, and if not offer some error
2075 # correction
2076 if self._kw_tol == 0:
2077 similars = difflib.get_close_matches(
2078 attr,
2079 self.cell._options.keys() + self.param._options.keys())
2080 if similars:
2081 raise UserWarning('Option "%s" not known! You mean "%s"?' %
2082 (attr, similars[0]))
2083 else:
2084 raise UserWarning('Option "%s" is not known!' % attr)
2085 else:
2086 warnings.warn('Option "%s" is not known - please set any new'
2087 ' options directly in the .cell or .param '
2088 'objects' % attr)
2089 return
2091 # here we know it must go into one of the component param or cell
2092 # so we first determine which one
2093 if attr in self.param._options.keys():
2094 comp = 'param'
2095 elif attr in self.cell._options.keys():
2096 comp = 'cell'
2097 else:
2098 raise UserWarning('Programming error: could not attach '
2099 'the keyword to an input file')
2101 self.__dict__[comp].__setattr__(attr, value)
2103 def merge_param(self, param, overwrite=True, ignore_internal_keys=False):
2104 """Parse a param file and merge it into the current parameters."""
2105 if isinstance(param, CastepParam):
2106 for key, option in param._options.items():
2107 if option.value is not None:
2108 self.param.__setattr__(key, option.value)
2109 return
2111 elif isinstance(param, str):
2112 param_file = open(param, 'r')
2113 _close = True
2115 else:
2116 # in this case we assume that we have a fileobj already, but check
2117 # for attributes in order to avoid extended EAFP blocks.
2118 param_file = param
2120 # look before you leap...
2121 attributes = ['name',
2122 'close'
2123 'readlines']
2125 for attr in attributes:
2126 if not hasattr(param_file, attr):
2127 raise TypeError('"param" is neither CastepParam nor str '
2128 'nor valid fileobj')
2130 param = param_file.name
2131 _close = False
2133 self, int_opts = read_param(fd=param_file, calc=self,
2134 get_interface_options=True)
2136 # Add the interface options
2137 for k, val in int_opts.items():
2138 if (k in self.internal_keys and not ignore_internal_keys):
2139 if val in _tf_table:
2140 val = _tf_table[val]
2141 self._opt[k] = val
2143 if _close:
2144 param_file.close()
2146 def dryrun_ok(self, dryrun_flag='-dryrun'):
2147 """Starts a CASTEP run with the -dryrun flag [default]
2148 in a temporary and check wether all variables are initialized
2149 correctly. This is recommended for every bigger simulation.
2150 """
2151 from ase.io.castep import write_param
2153 temp_dir = tempfile.mkdtemp()
2154 self._fetch_pspots(temp_dir)
2155 seed = 'dryrun'
2157 self._write_cell(os.path.join(temp_dir, '%s.cell' % seed),
2158 self.atoms, castep_cell=self.cell)
2159 # This part needs to be modified now that we rely on the new formats.py
2160 # interface
2161 if not os.path.isfile(os.path.join(temp_dir, '%s.cell' % seed)):
2162 warnings.warn('%s.cell not written - aborting dryrun' % seed)
2163 return
2164 write_param(os.path.join(temp_dir, '%s.param' % seed), self.param, )
2166 stdout, stderr = shell_stdouterr(('%s %s %s' % (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, '%s.castep' % seed))
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, '%s.0001.err' % seed)
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 os.environ.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('PSPs directory %s not found' % self._castep_pp_path)
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 'Warning: you have no PP specified for %s.' %
2243 species)
2244 if self.cell.species_pot.value:
2245 for (species, pspot) in pspots.items():
2246 orig_pspot_file = os.path.join(self._castep_pp_path, pspot)
2247 cp_pspot_file = os.path.join(directory, pspot)
2248 if (os.path.exists(orig_pspot_file)
2249 and not os.path.exists(cp_pspot_file)):
2250 if self._copy_pspots:
2251 shutil.copy(orig_pspot_file, directory)
2252 elif self._link_pspots:
2253 os.symlink(orig_pspot_file, cp_pspot_file)
2254 else:
2255 if self._pedantic:
2256 warnings.warn(ppwarning)
2259ppwarning = ('Warning: PP files have neither been '
2260 'linked nor copied to the working directory. Make '
2261 'sure to set the evironment variable PSPOT_DIR '
2262 'accordingly!')
2265def get_castep_version(castep_command):
2266 """This returns the version number as printed in the CASTEP banner.
2267 For newer CASTEP versions ( > 6.1) the --version command line option
2268 has been added; this will be attempted first.
2269 """
2270 import tempfile
2271 with tempfile.TemporaryDirectory() as temp_dir:
2272 return _get_castep_version(castep_command, temp_dir)
2275def _get_castep_version(castep_command, temp_dir):
2276 jname = 'dummy_jobname'
2277 stdout, stderr = '', ''
2278 fallback_version = 16. # CASTEP 16.0 and 16.1 report version wrongly
2279 try:
2280 stdout, stderr = subprocess.Popen(
2281 castep_command.split() + ['--version'],
2282 stderr=subprocess.PIPE,
2283 stdout=subprocess.PIPE, cwd=temp_dir,
2284 universal_newlines=True).communicate()
2285 if 'CASTEP version' not in stdout:
2286 stdout, stderr = subprocess.Popen(
2287 castep_command.split() + [jname],
2288 stderr=subprocess.PIPE,
2289 stdout=subprocess.PIPE, cwd=temp_dir,
2290 universal_newlines=True).communicate()
2291 except Exception: # XXX Which kind of exception?
2292 msg = ''
2293 msg += 'Could not determine the version of your CASTEP binary \n'
2294 msg += 'This usually means one of the following \n'
2295 msg += ' * you do not have CASTEP installed \n'
2296 msg += ' * you have not set the CASTEP_COMMAND to call it \n'
2297 msg += ' * you have provided a wrong CASTEP_COMMAND. \n'
2298 msg += ' Make sure it is in your PATH\n\n'
2299 msg += stdout
2300 msg += stderr
2301 raise CastepVersionError(msg)
2302 if 'CASTEP version' in stdout:
2303 output_txt = stdout.split('\n')
2304 version_re = re.compile(r'CASTEP version:\s*([0-9\.]*)')
2305 else:
2306 output = open(os.path.join(temp_dir, '%s.castep' % jname))
2307 output_txt = output.readlines()
2308 output.close()
2309 version_re = re.compile(r'(?<=CASTEP version )[0-9.]*')
2310 # shutil.rmtree(temp_dir)
2311 for line in output_txt:
2312 if 'CASTEP version' in line:
2313 try:
2314 return float(version_re.findall(line)[0])
2315 except ValueError:
2316 # Fallback for buggy --version on CASTEP 16.0, 16.1
2317 return fallback_version
2320def create_castep_keywords(castep_command, filename='castep_keywords.json',
2321 force_write=True, path='.', fetch_only=None):
2322 """This function allows to fetch all available keywords from stdout
2323 of an installed castep binary. It furthermore collects the documentation
2324 to harness the power of (ipython) inspection and type for some basic
2325 type checking of input. All information is stored in a JSON file that is
2326 not distributed by default to avoid breaking the license of CASTEP.
2327 """
2328 # Takes a while ...
2329 # Fetch all allowed parameters
2330 # fetch_only : only fetch that many parameters (for testsuite only)
2331 suffixes = ['cell', 'param']
2333 filepath = os.path.join(path, filename)
2335 if os.path.exists(filepath) and not force_write:
2336 warnings.warn('CASTEP Options Module file exists. '
2337 'You can overwrite it by calling '
2338 'python castep.py -f [CASTEP_COMMAND].')
2339 return False
2341 # Not saving directly to file her to prevent half-generated files
2342 # which will cause problems on future runs
2344 castep_version = get_castep_version(castep_command)
2346 help_all, _ = shell_stdouterr('%s -help all' % castep_command)
2348 # Filter out proper keywords
2349 try:
2350 # The old pattern does not math properly as in CASTEP as of v8.0 there
2351 # are some keywords for the semi-empircal dispersion correction (SEDC)
2352 # which also include numbers.
2353 if castep_version < 7.0:
2354 pattern = r'((?<=^ )[A-Z_]{2,}|(?<=^)[A-Z_]{2,})'
2355 else:
2356 pattern = r'((?<=^ )[A-Z_\d]{2,}|(?<=^)[A-Z_\d]{2,})'
2358 raw_options = re.findall(pattern, help_all, re.MULTILINE)
2359 except Exception:
2360 warnings.warn('Problem parsing: %s' % help_all)
2361 raise
2363 types = set()
2364 levels = set()
2366 processed_n = 0
2367 to_process = len(raw_options[:fetch_only])
2369 processed_options = {sf: {} for sf in suffixes}
2371 for o_i, option in enumerate(raw_options[:fetch_only]):
2372 doc, _ = shell_stdouterr('%s -help %s' % (castep_command, option))
2374 # Stand Back! I know regular expressions (http://xkcd.com/208/) :-)
2375 match = re.match(r'(?P<before_type>.*)Type: (?P<type>.+?)\s+'
2376 + r'Level: (?P<level>[^ ]+)\n\s*\n'
2377 + r'(?P<doc>.*?)(\n\s*\n|$)', doc, re.DOTALL)
2379 processed_n += 1
2381 if match is not None:
2382 match = match.groupdict()
2384 # JM: uncomment lines in following block to debug issues
2385 # with keyword assignment during extraction process from CASTEP
2386 suffix = None
2387 if re.findall(r'PARAMETERS keywords:\n\n\s?None found', doc):
2388 suffix = 'cell'
2389 if re.findall(r'CELL keywords:\n\n\s?None found', doc):
2390 suffix = 'param'
2391 if suffix is None:
2392 warnings.warn('%s -> not assigned to either'
2393 ' CELL or PARAMETERS keywords' % option)
2395 option = option.lower()
2396 mtyp = match.get('type', None)
2397 mlvl = match.get('level', None)
2398 mdoc = match.get('doc', None)
2400 if mtyp is None:
2401 warnings.warn('Found no type for %s' % option)
2402 continue
2403 if mlvl is None:
2404 warnings.warn('Found no level for %s' % option)
2405 continue
2406 if mdoc is None:
2407 warnings.warn('Found no doc string for %s' % option)
2408 continue
2410 types = types.union([mtyp])
2411 levels = levels.union([mlvl])
2413 processed_options[suffix][option] = {
2414 'keyword': option,
2415 'option_type': mtyp,
2416 'level': mlvl,
2417 'docstring': mdoc
2418 }
2420 processed_n += 1
2422 frac = (o_i + 1.0) / to_process
2423 sys.stdout.write('\rProcessed: [{0}] {1:>3.0f}%'.format(
2424 '#' * int(frac * 20) + ' '
2425 * (20 - int(frac * 20)),
2426 100 * frac))
2427 sys.stdout.flush()
2429 else:
2430 warnings.warn('create_castep_keywords: Could not process %s'
2431 % option)
2433 sys.stdout.write('\n')
2434 sys.stdout.flush()
2436 processed_options['types'] = list(types)
2437 processed_options['levels'] = list(levels)
2438 processed_options['castep_version'] = castep_version
2440 json.dump(processed_options, open(filepath, 'w'), indent=4)
2442 warnings.warn('CASTEP v%s, fetched %s keywords' %
2443 (castep_version, processed_n))
2444 return True
2447class CastepOption:
2448 """"A CASTEP option. It handles basic conversions from string to its value
2449 type."""
2451 default_convert_types = {
2452 'boolean (logical)': 'bool',
2453 'defined': 'bool',
2454 'string': 'str',
2455 'integer': 'int',
2456 'real': 'float',
2457 'integer vector': 'int_vector',
2458 'real vector': 'float_vector',
2459 'physical': 'float_physical',
2460 'block': 'block'
2461 }
2463 def __init__(self, keyword, level, option_type, value=None,
2464 docstring='No information available'):
2465 self.keyword = keyword
2466 self.level = level
2467 self.type = option_type
2468 self._value = value
2469 self.__doc__ = docstring
2471 @property
2472 def value(self):
2474 if self._value is not None:
2475 if self.type.lower() in ('integer vector', 'real vector',
2476 'physical'):
2477 return ' '.join(map(str, self._value))
2478 elif self.type.lower() in ('boolean (logical)', 'defined'):
2479 return str(self._value).upper()
2480 else:
2481 return str(self._value)
2483 @property
2484 def raw_value(self):
2485 # The value, not converted to a string
2486 return self._value
2488 @value.setter # type: ignore
2489 def value(self, val):
2491 if val is None:
2492 self.clear()
2493 return
2495 ctype = self.default_convert_types.get(self.type.lower(), 'str')
2496 typeparse = '_parse_%s' % ctype
2497 try:
2498 self._value = getattr(self, typeparse)(val)
2499 except ValueError:
2500 raise ConversionError(ctype, self.keyword, val)
2502 def clear(self):
2503 """Reset the value of the option to None again"""
2504 self._value = None
2506 @staticmethod
2507 def _parse_bool(value):
2508 try:
2509 value = _tf_table[str(value).strip().title()]
2510 except (KeyError, ValueError):
2511 raise ValueError()
2512 return value
2514 @staticmethod
2515 def _parse_str(value):
2516 value = str(value)
2517 return value
2519 @staticmethod
2520 def _parse_int(value):
2521 value = int(value)
2522 return value
2524 @staticmethod
2525 def _parse_float(value):
2526 value = float(value)
2527 return value
2529 @staticmethod
2530 def _parse_int_vector(value):
2531 # Accepts either a string or an actual list/numpy array of ints
2532 if isinstance(value, str):
2533 if ',' in value:
2534 value = value.replace(',', ' ')
2535 value = list(map(int, value.split()))
2537 value = np.array(value)
2539 if value.shape != (3,) or value.dtype != int:
2540 raise ValueError()
2542 return list(value)
2544 @staticmethod
2545 def _parse_float_vector(value):
2546 # Accepts either a string or an actual list/numpy array of floats
2547 if isinstance(value, str):
2548 if ',' in value:
2549 value = value.replace(',', ' ')
2550 value = list(map(float, value.split()))
2552 value = np.array(value) * 1.0
2554 if value.shape != (3,) or value.dtype != float:
2555 raise ValueError()
2557 return list(value)
2559 @staticmethod
2560 def _parse_float_physical(value):
2561 # If this is a string containing units, saves them
2562 if isinstance(value, str):
2563 value = value.split()
2565 try:
2566 l = len(value)
2567 except TypeError:
2568 l = 1
2569 value = [value]
2571 if l == 1:
2572 try:
2573 value = (float(value[0]), '')
2574 except (TypeError, ValueError):
2575 raise ValueError()
2576 elif l == 2:
2577 try:
2578 value = (float(value[0]), value[1])
2579 except (TypeError, ValueError, IndexError):
2580 raise ValueError()
2581 else:
2582 raise ValueError()
2584 return value
2586 @staticmethod
2587 def _parse_block(value):
2589 if isinstance(value, str):
2590 return value
2591 elif hasattr(value, '__getitem__'):
2592 return '\n'.join(value) # Arrays of lines
2593 else:
2594 raise ValueError()
2596 def __repr__(self):
2597 if self._value:
2598 expr = ('Option: {keyword}({type}, {level}):\n{_value}\n'
2599 ).format(**self.__dict__)
2600 else:
2601 expr = ('Option: {keyword}[unset]({type}, {level})'
2602 ).format(**self.__dict__)
2603 return expr
2605 def __eq__(self, other):
2606 if not isinstance(other, CastepOption):
2607 return False
2608 else:
2609 return self.__dict__ == other.__dict__
2612class CastepOptionDict:
2613 """A dictionary-like object to hold a set of options for .cell or .param
2614 files loaded from a dictionary, for the sake of validation.
2616 Replaces the old CastepCellDict and CastepParamDict that were defined in
2617 the castep_keywords.py file.
2618 """
2620 def __init__(self, options=None):
2621 object.__init__(self)
2622 self._options = {} # ComparableDict is not needed any more as
2623 # CastepOptions can be compared directly now
2624 for kw in options:
2625 opt = CastepOption(**options[kw])
2626 self._options[opt.keyword] = opt
2627 self.__dict__[opt.keyword] = opt
2630class CastepInputFile:
2632 """Master class for CastepParam and CastepCell to inherit from"""
2634 _keyword_conflicts: List[Set[str]] = []
2636 def __init__(self, options_dict=None, keyword_tolerance=1):
2637 object.__init__(self)
2639 if options_dict is None:
2640 options_dict = CastepOptionDict({})
2642 self._options = options_dict._options
2643 self.__dict__.update(self._options)
2644 # keyword_tolerance means how strict the checks on new attributes are
2645 # 0 = no new attributes allowed
2646 # 1 = new attributes allowed, warning given
2647 # 2 = new attributes allowed, silent
2648 self._perm = np.clip(keyword_tolerance, 0, 2)
2650 # Compile a dictionary for quick check of conflict sets
2651 self._conflict_dict = {
2652 kw: set(cset).difference({kw})
2653 for cset in self._keyword_conflicts for kw in cset}
2655 def __repr__(self):
2656 expr = ''
2657 is_default = True
2658 for key, option in sorted(self._options.items()):
2659 if option.value is not None:
2660 is_default = False
2661 expr += ('%20s : %s\n' % (key, option.value))
2662 if is_default:
2663 expr = 'Default\n'
2665 expr += 'Keyword tolerance: {0}'.format(self._perm)
2666 return expr
2668 def __setattr__(self, attr, value):
2670 # Hidden attributes are treated normally
2671 if attr.startswith('_'):
2672 self.__dict__[attr] = value
2673 return
2675 if attr not in self._options.keys():
2677 if self._perm > 0:
2678 # Do we consider it a string or a block?
2679 is_str = isinstance(value, str)
2680 is_block = False
2681 if ((hasattr(value, '__getitem__') and not is_str)
2682 or (is_str and len(value.split('\n')) > 1)):
2683 is_block = True
2685 if self._perm == 0:
2686 similars = difflib.get_close_matches(attr,
2687 self._options.keys())
2688 if similars:
2689 raise UserWarning(('Option "%s" not known! You mean "%s"?')
2690 % (attr, similars[0]))
2691 else:
2692 raise UserWarning('Option "%s" is not known!' % attr)
2693 elif self._perm == 1:
2694 warnings.warn(('Option "%s" is not known and will '
2695 'be added as a %s') % (attr,
2696 ('block' if is_block else
2697 'string')))
2698 attr = attr.lower()
2699 opt = CastepOption(keyword=attr, level='Unknown',
2700 option_type='block' if is_block else 'string')
2701 self._options[attr] = opt
2702 self.__dict__[attr] = opt
2703 else:
2704 attr = attr.lower()
2705 opt = self._options[attr]
2707 if not opt.type.lower() == 'block' and isinstance(value, str):
2708 value = value.replace(':', ' ')
2710 # If it is, use the appropriate parser, unless a custom one is defined
2711 attrparse = '_parse_%s' % attr.lower()
2713 # Check for any conflicts if the value is not None
2714 if not (value is None):
2715 cset = self._conflict_dict.get(attr.lower(), {})
2716 for c in cset:
2717 if (c in self._options and self._options[c].value):
2718 warnings.warn(
2719 'option "{attr}" conflicts with "{conflict}" in '
2720 'calculator. Setting "{conflict}" to '
2721 'None.'.format(attr=attr, conflict=c))
2722 self._options[c].value = None
2724 if hasattr(self, attrparse):
2725 self._options[attr].value = self.__getattribute__(attrparse)(value)
2726 else:
2727 self._options[attr].value = value
2729 def __getattr__(self, name):
2730 if name[0] == '_' or self._perm == 0:
2731 raise AttributeError()
2733 if self._perm == 1:
2734 warnings.warn('Option %s is not known, returning None' % (name))
2736 return CastepOption(keyword='none', level='Unknown',
2737 option_type='string', value=None)
2739 def get_attr_dict(self, raw=False, types=False):
2740 """Settings that go into .param file in a traditional dict"""
2742 attrdict = {k: o.raw_value if raw else o.value
2743 for k, o in self._options.items() if o.value is not None}
2745 if types:
2746 for key, val in attrdict.items():
2747 attrdict[key] = (val, self._options[key].type)
2749 return attrdict
2752class CastepParam(CastepInputFile):
2753 """CastepParam abstracts the settings that go into the .param file"""
2755 _keyword_conflicts = [{'cut_off_energy', 'basis_precision'}, ]
2757 def __init__(self, castep_keywords, keyword_tolerance=1):
2758 self._castep_version = castep_keywords.castep_version
2759 CastepInputFile.__init__(self, castep_keywords.CastepParamDict(),
2760 keyword_tolerance)
2762 @property
2763 def castep_version(self):
2764 return self._castep_version
2766 # .param specific parsers
2767 def _parse_reuse(self, value):
2768 if value is None:
2769 return None # Reset the value
2770 try:
2771 if self._options['continuation'].value:
2772 warnings.warn('Cannot set reuse if continuation is set, and '
2773 'vice versa. Set the other to None, if you want '
2774 'this setting.')
2775 return None
2776 except KeyError:
2777 pass
2778 return 'default' if (value is True) else str(value)
2780 def _parse_continuation(self, value):
2781 if value is None:
2782 return None # Reset the value
2783 try:
2784 if self._options['reuse'].value:
2785 warnings.warn('Cannot set reuse if continuation is set, and '
2786 'vice versa. Set the other to None, if you want '
2787 'this setting.')
2788 return None
2789 except KeyError:
2790 pass
2791 return 'default' if (value is True) else str(value)
2794class CastepCell(CastepInputFile):
2796 """CastepCell abstracts all setting that go into the .cell file"""
2798 _keyword_conflicts = [
2799 {'kpoint_mp_grid', 'kpoint_mp_spacing', 'kpoint_list',
2800 'kpoints_mp_grid', 'kpoints_mp_spacing', 'kpoints_list'},
2801 {'bs_kpoint_mp_grid',
2802 'bs_kpoint_mp_spacing',
2803 'bs_kpoint_list',
2804 'bs_kpoint_path',
2805 'bs_kpoints_mp_grid',
2806 'bs_kpoints_mp_spacing',
2807 'bs_kpoints_list',
2808 'bs_kpoints_path'},
2809 {'spectral_kpoint_mp_grid',
2810 'spectral_kpoint_mp_spacing',
2811 'spectral_kpoint_list',
2812 'spectral_kpoint_path',
2813 'spectral_kpoints_mp_grid',
2814 'spectral_kpoints_mp_spacing',
2815 'spectral_kpoints_list',
2816 'spectral_kpoints_path'},
2817 {'phonon_kpoint_mp_grid',
2818 'phonon_kpoint_mp_spacing',
2819 'phonon_kpoint_list',
2820 'phonon_kpoint_path',
2821 'phonon_kpoints_mp_grid',
2822 'phonon_kpoints_mp_spacing',
2823 'phonon_kpoints_list',
2824 'phonon_kpoints_path'},
2825 {'fine_phonon_kpoint_mp_grid',
2826 'fine_phonon_kpoint_mp_spacing',
2827 'fine_phonon_kpoint_list',
2828 'fine_phonon_kpoint_path'},
2829 {'magres_kpoint_mp_grid',
2830 'magres_kpoint_mp_spacing',
2831 'magres_kpoint_list',
2832 'magres_kpoint_path'},
2833 {'elnes_kpoint_mp_grid',
2834 'elnes_kpoint_mp_spacing',
2835 'elnes_kpoint_list',
2836 'elnes_kpoint_path'},
2837 {'optics_kpoint_mp_grid',
2838 'optics_kpoint_mp_spacing',
2839 'optics_kpoint_list',
2840 'optics_kpoint_path'},
2841 {'supercell_kpoint_mp_grid',
2842 'supercell_kpoint_mp_spacing',
2843 'supercell_kpoint_list',
2844 'supercell_kpoint_path'}, ]
2846 def __init__(self, castep_keywords, keyword_tolerance=1):
2847 self._castep_version = castep_keywords.castep_version
2848 CastepInputFile.__init__(self, castep_keywords.CastepCellDict(),
2849 keyword_tolerance)
2851 @property
2852 def castep_version(self):
2853 return self._castep_version
2855 # .cell specific parsers
2856 def _parse_species_pot(self, value):
2858 # Single tuple
2859 if isinstance(value, tuple) and len(value) == 2:
2860 value = [value]
2861 # List of tuples
2862 if hasattr(value, '__getitem__'):
2863 pspots = [tuple(map(str.strip, x)) for x in value]
2864 if not all(map(lambda x: len(x) == 2, value)):
2865 warnings.warn(
2866 'Please specify pseudopotentials in python as '
2867 'a tuple or a list of tuples formatted like: '
2868 '(species, file), e.g. ("O", "path-to/O_OTFG.usp") '
2869 'Anything else will be ignored')
2870 return None
2872 text_block = self._options['species_pot'].value
2874 text_block = text_block if text_block else ''
2875 # Remove any duplicates
2876 for pp in pspots:
2877 text_block = re.sub(r'\n?\s*%s\s+.*' % pp[0], '', text_block)
2878 if pp[1]:
2879 text_block += '\n%s %s' % pp
2881 return text_block
2883 def _parse_symmetry_ops(self, value):
2884 if not isinstance(value, tuple) \
2885 or not len(value) == 2 \
2886 or not value[0].shape[1:] == (3, 3) \
2887 or not value[1].shape[1:] == (3,) \
2888 or not value[0].shape[0] == value[1].shape[0]:
2889 warnings.warn('Invalid symmetry_ops block, skipping')
2890 return
2891 # Now on to print...
2892 text_block = ''
2893 for op_i, (op_rot, op_tranls) in enumerate(zip(*value)):
2894 text_block += '\n'.join([' '.join([str(x) for x in row])
2895 for row in op_rot])
2896 text_block += '\n'
2897 text_block += ' '.join([str(x) for x in op_tranls])
2898 text_block += '\n\n'
2900 return text_block
2902 def _parse_positions_abs_intermediate(self, value):
2903 return _parse_tss_block(value)
2905 def _parse_positions_abs_product(self, value):
2906 return _parse_tss_block(value)
2908 def _parse_positions_frac_intermediate(self, value):
2909 return _parse_tss_block(value, True)
2911 def _parse_positions_frac_product(self, value):
2912 return _parse_tss_block(value, True)
2915CastepKeywords = namedtuple('CastepKeywords',
2916 ['CastepParamDict', 'CastepCellDict',
2917 'types', 'levels', 'castep_version'])
2919# We keep this just for naming consistency with older versions
2922def make_cell_dict(data=None):
2924 data = data if data is not None else {}
2926 class CastepCellDict(CastepOptionDict):
2927 def __init__(self):
2928 CastepOptionDict.__init__(self, data)
2930 return CastepCellDict
2933def make_param_dict(data=None):
2935 data = data if data is not None else {}
2937 class CastepParamDict(CastepOptionDict):
2938 def __init__(self):
2939 CastepOptionDict.__init__(self, data)
2941 return CastepParamDict
2944class CastepVersionError(Exception):
2945 """No special behaviour, works to signal when Castep can not be found"""
2946 pass
2949class ConversionError(Exception):
2951 """Print customized error for options that are not converted correctly
2952 and point out that they are maybe not implemented, yet"""
2954 def __init__(self, key_type, attr, value):
2955 Exception.__init__(self)
2956 self.key_type = key_type
2957 self.value = value
2958 self.attr = attr
2960 def __str__(self):
2961 return 'Could not convert %s = %s to %s\n' \
2962 % (self.attr, self.value, self.key_type) \
2963 + 'This means you either tried to set a value of the wrong\n'\
2964 + 'type or this keyword needs some special care. Please feel\n'\
2965 + 'to add it to the corresponding __setattr__ method and send\n'\
2966 + 'the patch to %s, so we can all benefit.' % (contact_email)
2969def get_castep_pp_path(castep_pp_path=''):
2970 """Abstract the quest for a CASTEP PSP directory."""
2971 if castep_pp_path:
2972 return os.path.abspath(os.path.expanduser(castep_pp_path))
2973 elif 'PSPOT_DIR' in os.environ:
2974 return os.environ['PSPOT_DIR']
2975 elif 'CASTEP_PP_PATH' in os.environ:
2976 return os.environ['CASTEP_PP_PATH']
2977 else:
2978 return os.path.abspath('.')
2981def get_castep_command(castep_command=''):
2982 """Abstract the quest for a castep_command string."""
2983 if castep_command:
2984 return castep_command
2985 elif 'CASTEP_COMMAND' in os.environ:
2986 return os.environ['CASTEP_COMMAND']
2987 else:
2988 return 'castep'
2991def shell_stdouterr(raw_command, cwd=None):
2992 """Abstracts the standard call of the commandline, when
2993 we are only interested in the stdout and stderr
2994 """
2995 stdout, stderr = subprocess.Popen(raw_command,
2996 stdout=subprocess.PIPE,
2997 stderr=subprocess.PIPE,
2998 universal_newlines=True,
2999 shell=True, cwd=cwd).communicate()
3000 return stdout.strip(), stderr.strip()
3003def import_castep_keywords(castep_command='',
3004 filename='castep_keywords.json',
3005 path='.'):
3006 """Search for castep keywords JSON in multiple paths"""
3008 config_paths = ('~/.ase', '~/.config/ase')
3009 searchpaths = [path] + [os.path.expanduser(config_path)
3010 for config_path in config_paths]
3011 try:
3012 keywords_file = sum([glob.glob(os.path.join(sp, filename))
3013 for sp in searchpaths], [])[0]
3014 except IndexError:
3015 warnings.warn("""Generating CASTEP keywords JSON file... hang on.
3016 The CASTEP keywords JSON file contains abstractions for CASTEP input
3017 parameters (for both .cell and .param input files), including some
3018 format checks and descriptions. The latter are extracted from the
3019 internal online help facility of a CASTEP binary, thus allowing to
3020 easily keep the calculator synchronized with (different versions of)
3021 the CASTEP code. Consequently, avoiding licensing issues (CASTEP is
3022 distributed commercially by Biovia), we consider it wise not to
3023 provide the file in the first place.""")
3024 create_castep_keywords(get_castep_command(castep_command),
3025 filename=filename, path=path)
3026 keywords_file = Path(path).absolute() / filename
3028 warnings.warn(
3029 f'Stored castep keywords dictionary as {keywords_file}. '
3030 f'Copy it to {Path(config_paths[0]).expanduser() / filename} for '
3031 r'user installation.')
3033 # Now create the castep_keywords object proper
3034 with open(keywords_file) as fd:
3035 kwdata = json.load(fd)
3037 # This is a bit awkward, but it's necessary for backwards compatibility
3038 param_dict = make_param_dict(kwdata['param'])
3039 cell_dict = make_cell_dict(kwdata['cell'])
3041 castep_keywords = CastepKeywords(param_dict, cell_dict,
3042 kwdata['types'], kwdata['levels'],
3043 kwdata['castep_version'])
3045 return castep_keywords
3048if __name__ == '__main__':
3049 warnings.warn(
3050 'When called directly this calculator will fetch all available '
3051 'keywords from the binarys help function into a '
3052 'castep_keywords.json in the current directory %s '
3053 'For system wide usage, it can be copied into an ase installation '
3054 'at ASE/calculators. '
3055 'This castep_keywords.json usually only needs to be generated once '
3056 'for a CASTEP binary/CASTEP version.' % os.getcwd())
3058 import optparse
3059 parser = optparse.OptionParser()
3060 parser.add_option(
3061 '-f', '--force-write', dest='force_write',
3062 help='Force overwriting existing castep_keywords.json', default=False,
3063 action='store_true')
3064 (options, args) = parser.parse_args()
3066 if args:
3067 opt_castep_command = ''.join(args)
3068 else:
3069 opt_castep_command = ''
3070 generated = create_castep_keywords(get_castep_command(opt_castep_command),
3071 force_write=options.force_write)
3073 if generated:
3074 try:
3075 with open('castep_keywords.json') as fd:
3076 json.load(fd)
3077 except Exception as e:
3078 warnings.warn(
3079 '%s Ooops, something went wrong with the CASTEP keywords' % e)
3080 else:
3081 warnings.warn('Import works. Looking good!')