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