Coverage for /builds/debichem-team/python-ase/ase/calculators/calculator.py: 83.85%
545 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-03-06 04:00 +0000
1import copy
2import os
3import shlex
4import subprocess
5import warnings
6from abc import abstractmethod
7from dataclasses import dataclass, field
8from math import pi, sqrt
9from pathlib import Path
10from typing import Any, Dict, List, Optional, Sequence, Set, Union
12import numpy as np
14from ase.calculators.abc import GetPropertiesMixin
15from ase.cell import Cell
16from ase.config import cfg as _cfg
17from ase.outputs import Properties, all_outputs
18from ase.utils import deprecated, jsonable
20from .names import names
23class CalculatorError(RuntimeError):
24 """Base class of error types related to ASE calculators."""
27class CalculatorSetupError(CalculatorError):
28 """Calculation cannot be performed with the given parameters.
30 Reasons to raise this error are:
31 * The calculator is not properly configured
32 (missing executable, environment variables, ...)
33 * The given atoms object is not supported
34 * Calculator parameters are unsupported
36 Typically raised before a calculation."""
39class EnvironmentError(CalculatorSetupError):
40 """Raised if calculator is not properly set up with ASE.
42 May be missing an executable or environment variables."""
45class InputError(CalculatorSetupError):
46 """Raised if inputs given to the calculator were incorrect.
48 Bad input keywords or values, or missing pseudopotentials.
50 This may be raised before or during calculation, depending on
51 when the problem is detected."""
54class CalculationFailed(CalculatorError):
55 """Calculation failed unexpectedly.
57 Reasons to raise this error are:
58 * Calculation did not converge
59 * Calculation ran out of memory
60 * Segmentation fault or other abnormal termination
61 * Arithmetic trouble (singular matrices, NaN, ...)
63 Typically raised during calculation."""
66class SCFError(CalculationFailed):
67 """SCF loop did not converge."""
70class ReadError(CalculatorError):
71 """Unexpected irrecoverable error while reading calculation results."""
74class PropertyNotImplementedError(NotImplementedError):
75 """Raised if a calculator does not implement the requested property."""
78class PropertyNotPresent(CalculatorError):
79 """Requested property is missing.
81 Maybe it was never calculated, or for some reason was not extracted
82 with the rest of the results, without being a fatal ReadError."""
85def compare_atoms(atoms1, atoms2, tol=1e-15, excluded_properties=None):
86 """Check for system changes since last calculation. Properties in
87 ``excluded_properties`` are not checked."""
88 if atoms1 is None:
89 system_changes = all_changes[:]
90 else:
91 system_changes = []
93 properties_to_check = set(all_changes)
94 if excluded_properties:
95 properties_to_check -= set(excluded_properties)
97 # Check properties that aren't in Atoms.arrays but are attributes of
98 # Atoms objects
99 for prop in ['cell', 'pbc']:
100 if prop in properties_to_check:
101 properties_to_check.remove(prop)
102 if not equal(
103 getattr(atoms1, prop), getattr(atoms2, prop), atol=tol
104 ):
105 system_changes.append(prop)
107 arrays1 = set(atoms1.arrays)
108 arrays2 = set(atoms2.arrays)
110 # Add any properties that are only in atoms1.arrays or only in
111 # atoms2.arrays (and aren't excluded). Note that if, e.g. arrays1 has
112 # `initial_charges` which is merely zeros and arrays2 does not have
113 # this array, we'll still assume that the system has changed. However,
114 # this should only occur rarely.
115 system_changes += properties_to_check & (arrays1 ^ arrays2)
117 # Finally, check all of the non-excluded properties shared by the atoms
118 # arrays.
119 for prop in properties_to_check & arrays1 & arrays2:
120 if not equal(atoms1.arrays[prop], atoms2.arrays[prop], atol=tol):
121 system_changes.append(prop)
123 return system_changes
126all_properties = [
127 'energy',
128 'forces',
129 'stress',
130 'stresses',
131 'dipole',
132 'charges',
133 'magmom',
134 'magmoms',
135 'free_energy',
136 'energies',
137 'dielectric_tensor',
138 'born_effective_charges',
139 'polarization',
140]
143all_changes = [
144 'positions',
145 'numbers',
146 'cell',
147 'pbc',
148 'initial_charges',
149 'initial_magmoms',
150]
153special = {
154 'cp2k': 'CP2K',
155 'demonnano': 'DemonNano',
156 'dftd3': 'DFTD3',
157 'dmol': 'DMol3',
158 'eam': 'EAM',
159 'elk': 'ELK',
160 'emt': 'EMT',
161 'exciting': 'ExcitingGroundStateCalculator',
162 'crystal': 'CRYSTAL',
163 'ff': 'ForceField',
164 'gamess_us': 'GAMESSUS',
165 'gulp': 'GULP',
166 'kim': 'KIM',
167 'lammpsrun': 'LAMMPS',
168 'lammpslib': 'LAMMPSlib',
169 'lj': 'LennardJones',
170 'mopac': 'MOPAC',
171 'morse': 'MorsePotential',
172 'nwchem': 'NWChem',
173 'openmx': 'OpenMX',
174 'orca': 'ORCA',
175 'qchem': 'QChem',
176 'tip3p': 'TIP3P',
177 'tip4p': 'TIP4P',
178}
181external_calculators = {}
184def register_calculator_class(name, cls):
185 """Add the class into the database."""
186 assert name not in external_calculators
187 external_calculators[name] = cls
188 names.append(name)
189 names.sort()
192def get_calculator_class(name):
193 """Return calculator class."""
194 if name == 'asap':
195 from asap3 import EMT as Calculator
196 elif name == 'gpaw':
197 from gpaw import GPAW as Calculator
198 elif name == 'hotbit':
199 from hotbit import Calculator
200 elif name == 'vasp2':
201 from ase.calculators.vasp import Vasp2 as Calculator
202 elif name == 'ace':
203 from ase.calculators.acemolecule import ACE as Calculator
204 elif name == 'Psi4':
205 from ase.calculators.psi4 import Psi4 as Calculator
206 elif name in external_calculators:
207 Calculator = external_calculators[name]
208 else:
209 classname = special.get(name, name.title())
210 module = __import__('ase.calculators.' + name, {}, None, [classname])
211 Calculator = getattr(module, classname)
212 return Calculator
215def equal(a, b, tol=None, rtol=None, atol=None):
216 """ndarray-enabled comparison function."""
217 # XXX Known bugs:
218 # * Comparing cell objects (pbc not part of array representation)
219 # * Infinite recursion for cyclic dicts
220 # * Can of worms is open
221 if tol is not None:
222 msg = 'Use `equal(a, b, rtol=..., atol=...)` instead of `tol=...`'
223 warnings.warn(msg, DeprecationWarning)
224 assert (
225 rtol is None and atol is None
226 ), 'Do not use deprecated `tol` with `atol` and/or `rtol`'
227 rtol = tol
228 atol = tol
230 a_is_dict = isinstance(a, dict)
231 b_is_dict = isinstance(b, dict)
232 if a_is_dict or b_is_dict:
233 # Check that both a and b are dicts
234 if not (a_is_dict and b_is_dict):
235 return False
236 if a.keys() != b.keys():
237 return False
238 return all(equal(a[key], b[key], rtol=rtol, atol=atol) for key in a)
240 if np.shape(a) != np.shape(b):
241 return False
243 if rtol is None and atol is None:
244 return np.array_equal(a, b)
246 if rtol is None:
247 rtol = 0
248 if atol is None:
249 atol = 0
251 return np.allclose(a, b, rtol=rtol, atol=atol)
254def kptdensity2monkhorstpack(atoms, kptdensity=3.5, even=True):
255 """Convert k-point density to Monkhorst-Pack grid size.
257 atoms: Atoms object
258 Contains unit cell and information about boundary conditions.
259 kptdensity: float
260 Required k-point density. Default value is 3.5 point per Ang^-1.
261 even: bool
262 Round up to even numbers.
263 """
265 recipcell = atoms.cell.reciprocal()
266 kpts = []
267 for i in range(3):
268 if atoms.pbc[i]:
269 k = 2 * pi * sqrt((recipcell[i] ** 2).sum()) * kptdensity
270 if even:
271 kpts.append(2 * int(np.ceil(k / 2)))
272 else:
273 kpts.append(int(np.ceil(k)))
274 else:
275 kpts.append(1)
276 return np.array(kpts)
279def kpts2mp(atoms, kpts, even=False):
280 if kpts is None:
281 return np.array([1, 1, 1])
282 if isinstance(kpts, (float, int)):
283 return kptdensity2monkhorstpack(atoms, kpts, even)
284 else:
285 return kpts
288def kpts2sizeandoffsets(
289 size=None, density=None, gamma=None, even=None, atoms=None
290):
291 """Helper function for selecting k-points.
293 Use either size or density.
295 size: 3 ints
296 Number of k-points.
297 density: float
298 K-point density in units of k-points per Ang^-1.
299 gamma: None or bool
300 Should the Gamma-point be included? Yes / no / don't care:
301 True / False / None.
302 even: None or bool
303 Should the number of k-points be even? Yes / no / don't care:
304 True / False / None.
305 atoms: Atoms object
306 Needed for calculating k-point density.
308 """
310 if size is not None and density is not None:
311 raise ValueError(
312 'Cannot specify k-point mesh size and ' 'density simultaneously'
313 )
314 elif density is not None and atoms is None:
315 raise ValueError(
316 'Cannot set k-points from "density" unless '
317 'Atoms are provided (need BZ dimensions).'
318 )
320 if size is None:
321 if density is None:
322 size = [1, 1, 1]
323 else:
324 size = kptdensity2monkhorstpack(atoms, density, None)
326 # Not using the rounding from kptdensity2monkhorstpack as it doesn't do
327 # rounding to odd numbers
328 if even is not None:
329 size = np.array(size)
330 remainder = size % 2
331 if even:
332 size += remainder
333 else: # Round up to odd numbers
334 size += 1 - remainder
336 offsets = [0, 0, 0]
337 if atoms is None:
338 pbc = [True, True, True]
339 else:
340 pbc = atoms.pbc
342 if gamma is not None:
343 for i, s in enumerate(size):
344 if pbc[i] and s % 2 != bool(gamma):
345 offsets[i] = 0.5 / s
347 return size, offsets
350@jsonable('kpoints')
351class KPoints:
352 def __init__(self, kpts=None):
353 if kpts is None:
354 kpts = np.zeros((1, 3))
355 self.kpts = kpts
357 def todict(self):
358 return vars(self)
361def kpts2kpts(kpts, atoms=None):
362 from ase.dft.kpoints import monkhorst_pack
364 if kpts is None:
365 return KPoints()
367 if hasattr(kpts, 'kpts'):
368 return kpts
370 if isinstance(kpts, dict):
371 if 'kpts' in kpts:
372 return KPoints(kpts['kpts'])
373 if 'path' in kpts:
374 cell = Cell.ascell(atoms.cell)
375 return cell.bandpath(pbc=atoms.pbc, **kpts)
376 size, offsets = kpts2sizeandoffsets(atoms=atoms, **kpts)
377 return KPoints(monkhorst_pack(size) + offsets)
379 if isinstance(kpts[0], int):
380 return KPoints(monkhorst_pack(kpts))
382 return KPoints(np.array(kpts))
385def kpts2ndarray(kpts, atoms=None):
386 """Convert kpts keyword to 2-d ndarray of scaled k-points."""
387 return kpts2kpts(kpts, atoms=atoms).kpts
390class EigenvalOccupationMixin:
391 """Define 'eigenvalues' and 'occupations' properties on class.
393 eigenvalues and occupations will be arrays of shape (spin, kpts, nbands).
395 Classes must implement the old-fashioned get_eigenvalues and
396 get_occupations methods."""
398 # We should maybe deprecate this and rely on the new
399 # Properties object for eigenvalues/occupations.
401 @property
402 def eigenvalues(self):
403 return self._propwrapper().eigenvalues
405 @property
406 def occupations(self):
407 return self._propwrapper().occupations
409 def _propwrapper(self):
410 from ase.calculator.singlepoint import OutputPropertyWrapper
412 return OutputPropertyWrapper(self)
415class Parameters(dict):
416 """Dictionary for parameters.
418 Special feature: If param is a Parameters instance, then param.xc
419 is a shorthand for param['xc'].
420 """
422 def __getattr__(self, key):
423 if key not in self:
424 return dict.__getattribute__(self, key)
425 return self[key]
427 def __setattr__(self, key, value):
428 self[key] = value
430 @classmethod
431 def read(cls, filename):
432 """Read parameters from file."""
433 # We use ast to evaluate literals, avoiding eval()
434 # for security reasons.
435 import ast
437 with open(filename) as fd:
438 txt = fd.read().strip()
439 assert txt.startswith('dict(')
440 assert txt.endswith(')')
441 txt = txt[5:-1]
443 # The tostring() representation "dict(...)" is not actually
444 # a literal, so we manually parse that along with the other
445 # formatting that we did manually:
446 dct = {}
447 for line in txt.splitlines():
448 key, val = line.split('=', 1)
449 key = key.strip()
450 val = val.strip()
451 if val[-1] == ',':
452 val = val[:-1]
453 dct[key] = ast.literal_eval(val)
455 parameters = cls(dct)
456 return parameters
458 def tostring(self):
459 keys = sorted(self)
460 return (
461 'dict('
462 + ',\n '.join(f'{key}={self[key]!r}' for key in keys)
463 + ')\n'
464 )
466 def write(self, filename):
467 Path(filename).write_text(self.tostring())
470class BaseCalculator(GetPropertiesMixin):
471 implemented_properties: List[str] = []
472 'Properties calculator can handle (energy, forces, ...)'
474 # Placeholder object for deprecated arguments. Let deprecated keywords
475 # default to _deprecated and then issue a warning if the user passed
476 # any other object (such as None).
477 _deprecated = object()
479 def __init__(self, parameters=None, use_cache=True):
480 if parameters is None:
481 parameters = {}
483 self.parameters = dict(parameters)
484 self.atoms = None
485 self.results = {}
486 self.use_cache = use_cache
488 def calculate_properties(self, atoms, properties):
489 """This method is experimental; currently for internal use."""
490 for name in properties:
491 if name not in all_outputs:
492 raise ValueError(f'No such property: {name}')
494 # We ignore system changes for now.
495 self.calculate(atoms, properties, system_changes=all_changes)
497 props = self.export_properties()
499 for name in properties:
500 if name not in props:
501 raise PropertyNotPresent(name)
502 return props
504 @abstractmethod
505 def calculate(self, atoms, properties, system_changes):
506 ...
508 def check_state(self, atoms, tol=1e-15):
509 """Check for any system changes since last calculation."""
510 if self.use_cache:
511 return compare_atoms(self.atoms, atoms, tol=tol)
512 else:
513 return all_changes
515 def get_property(self, name, atoms=None, allow_calculation=True):
516 if name not in self.implemented_properties:
517 raise PropertyNotImplementedError(
518 f'{name} property not implemented'
519 )
521 if atoms is None:
522 atoms = self.atoms
523 system_changes = []
524 else:
525 system_changes = self.check_state(atoms)
527 if system_changes:
528 self.atoms = None
529 self.results = {}
531 if name not in self.results:
532 if not allow_calculation:
533 return None
535 if self.use_cache:
536 self.atoms = atoms.copy()
538 self.calculate(atoms, [name], system_changes)
540 if name not in self.results:
541 # For some reason the calculator was not able to do what we want,
542 # and that is OK.
543 raise PropertyNotImplementedError(
544 '{} not present in this ' 'calculation'.format(name)
545 )
547 result = self.results[name]
548 if isinstance(result, np.ndarray):
549 result = result.copy()
550 return result
552 def calculation_required(self, atoms, properties):
553 assert not isinstance(properties, str)
554 system_changes = self.check_state(atoms)
555 if system_changes:
556 return True
557 for name in properties:
558 if name not in self.results:
559 return True
560 return False
562 def export_properties(self):
563 return Properties(self.results)
565 def _get_name(self) -> str: # child class can override this
566 return self.__class__.__name__.lower()
568 @property
569 def name(self) -> str:
570 return self._get_name()
573class Calculator(BaseCalculator):
574 """Base-class for all ASE calculators.
576 A calculator must raise PropertyNotImplementedError if asked for a
577 property that it can't calculate. So, if calculation of the
578 stress tensor has not been implemented, get_stress(atoms) should
579 raise PropertyNotImplementedError. This can be achieved simply by not
580 including the string 'stress' in the list implemented_properties
581 which is a class member. These are the names of the standard
582 properties: 'energy', 'forces', 'stress', 'dipole', 'charges',
583 'magmom' and 'magmoms'.
584 """
586 default_parameters: Dict[str, Any] = {}
587 'Default parameters'
589 ignored_changes: Set[str] = set()
590 'Properties of Atoms which we ignore for the purposes of cache '
591 'invalidation with check_state().'
593 discard_results_on_any_change = False
594 'Whether we purge the results following any change in the set() method. '
595 'Most (file I/O) calculators will probably want this.'
597 def __init__(
598 self,
599 restart=None,
600 ignore_bad_restart_file=BaseCalculator._deprecated,
601 label=None,
602 atoms=None,
603 directory='.',
604 **kwargs,
605 ):
606 """Basic calculator implementation.
608 restart: str
609 Prefix for restart file. May contain a directory. Default
610 is None: don't restart.
611 ignore_bad_restart_file: bool
612 Deprecated, please do not use.
613 Passing more than one positional argument to Calculator()
614 is deprecated and will stop working in the future.
615 Ignore broken or missing restart file. By default, it is an
616 error if the restart file is missing or broken.
617 directory: str or PurePath
618 Working directory in which to read and write files and
619 perform calculations.
620 label: str
621 Name used for all files. Not supported by all calculators.
622 May contain a directory, but please use the directory parameter
623 for that instead.
624 atoms: Atoms object
625 Optional Atoms object to which the calculator will be
626 attached. When restarting, atoms will get its positions and
627 unit-cell updated from file.
628 """
629 self.atoms = None # copy of atoms object from last calculation
630 self.results = {} # calculated properties (energy, forces, ...)
631 self.parameters = None # calculational parameters
632 self._directory = None # Initialize
634 if ignore_bad_restart_file is self._deprecated:
635 ignore_bad_restart_file = False
636 else:
637 warnings.warn(
638 FutureWarning(
639 'The keyword "ignore_bad_restart_file" is deprecated and '
640 'will be removed in a future version of ASE. Passing more '
641 'than one positional argument to Calculator is also '
642 'deprecated and will stop functioning in the future. '
643 'Please pass arguments by keyword (key=value) except '
644 'optionally the "restart" keyword.'
645 )
646 )
648 if restart is not None:
649 try:
650 self.read(restart) # read parameters, atoms and results
651 except ReadError:
652 if ignore_bad_restart_file:
653 self.reset()
654 else:
655 raise
657 self.directory = directory
658 self.prefix = None
659 if label is not None:
660 if self.directory == '.' and '/' in label:
661 # We specified directory in label, and nothing in the diretory
662 # key
663 self.label = label
664 elif '/' not in label:
665 # We specified our directory in the directory keyword
666 # or not at all
667 self.label = '/'.join((self.directory, label))
668 else:
669 raise ValueError(
670 'Directory redundantly specified though '
671 'directory="{}" and label="{}". '
672 'Please omit "/" in label.'.format(self.directory, label)
673 )
675 if self.parameters is None:
676 # Use default parameters if they were not read from file:
677 self.parameters = self.get_default_parameters()
679 if atoms is not None:
680 atoms.calc = self
681 if self.atoms is not None:
682 # Atoms were read from file. Update atoms:
683 if not (
684 equal(atoms.numbers, self.atoms.numbers)
685 and (atoms.pbc == self.atoms.pbc).all()
686 ):
687 raise CalculatorError('Atoms not compatible with file')
688 atoms.positions = self.atoms.positions
689 atoms.cell = self.atoms.cell
691 self.set(**kwargs)
693 if not hasattr(self, 'get_spin_polarized'):
694 self.get_spin_polarized = self._deprecated_get_spin_polarized
695 # XXX We are very naughty and do not call super constructor!
697 # For historical reasons we have a particular caching protocol.
698 # We disable the superclass' optional cache.
699 self.use_cache = False
701 @property
702 def directory(self) -> str:
703 return self._directory
705 @directory.setter
706 def directory(self, directory: Union[str, os.PathLike]):
707 self._directory = str(Path(directory)) # Normalize path.
709 @property
710 def label(self):
711 if self.directory == '.':
712 return self.prefix
714 # Generally, label ~ directory/prefix
715 #
716 # We use '/' rather than os.pathsep because
717 # 1) directory/prefix does not represent any actual path
718 # 2) We want the same string to work the same on all platforms
719 if self.prefix is None:
720 return self.directory + '/'
722 return f'{self.directory}/{self.prefix}'
724 @label.setter
725 def label(self, label):
726 if label is None:
727 self.directory = '.'
728 self.prefix = None
729 return
731 tokens = label.rsplit('/', 1)
732 if len(tokens) == 2:
733 directory, prefix = tokens
734 else:
735 assert len(tokens) == 1
736 directory = '.'
737 prefix = tokens[0]
738 if prefix == '':
739 prefix = None
740 self.directory = directory
741 self.prefix = prefix
743 def set_label(self, label):
744 """Set label and convert label to directory and prefix.
746 Examples:
748 * label='abc': (directory='.', prefix='abc')
749 * label='dir1/abc': (directory='dir1', prefix='abc')
750 * label=None: (directory='.', prefix=None)
751 """
752 self.label = label
754 def get_default_parameters(self):
755 return Parameters(copy.deepcopy(self.default_parameters))
757 def todict(self, skip_default=True):
758 defaults = self.get_default_parameters()
759 dct = {}
760 for key, value in self.parameters.items():
761 if hasattr(value, 'todict'):
762 value = value.todict()
763 if skip_default:
764 default = defaults.get(key, '_no_default_')
765 if default != '_no_default_' and equal(value, default):
766 continue
767 dct[key] = value
768 return dct
770 def reset(self):
771 """Clear all information from old calculation."""
773 self.atoms = None
774 self.results = {}
776 def read(self, label):
777 """Read atoms, parameters and calculated properties from output file.
779 Read result from self.label file. Raise ReadError if the file
780 is not there. If the file is corrupted or contains an error
781 message from the calculation, a ReadError should also be
782 raised. In case of succes, these attributes must set:
784 atoms: Atoms object
785 The state of the atoms from last calculation.
786 parameters: Parameters object
787 The parameter dictionary.
788 results: dict
789 Calculated properties like energy and forces.
791 The FileIOCalculator.read() method will typically read atoms
792 and parameters and get the results dict by calling the
793 read_results() method."""
795 self.set_label(label)
797 def get_atoms(self):
798 if self.atoms is None:
799 raise ValueError('Calculator has no atoms')
800 atoms = self.atoms.copy()
801 atoms.calc = self
802 return atoms
804 @classmethod
805 def read_atoms(cls, restart, **kwargs):
806 return cls(restart=restart, label=restart, **kwargs).get_atoms()
808 def set(self, **kwargs):
809 """Set parameters like set(key1=value1, key2=value2, ...).
811 A dictionary containing the parameters that have been changed
812 is returned.
814 Subclasses must implement a set() method that will look at the
815 chaneged parameters and decide if a call to reset() is needed.
816 If the changed parameters are harmless, like a change in
817 verbosity, then there is no need to call reset().
819 The special keyword 'parameters' can be used to read
820 parameters from a file."""
822 if 'parameters' in kwargs:
823 filename = kwargs.pop('parameters')
824 parameters = Parameters.read(filename)
825 parameters.update(kwargs)
826 kwargs = parameters
828 changed_parameters = {}
830 for key, value in kwargs.items():
831 oldvalue = self.parameters.get(key)
832 if key not in self.parameters or not equal(value, oldvalue):
833 changed_parameters[key] = value
834 self.parameters[key] = value
836 if self.discard_results_on_any_change and changed_parameters:
837 self.reset()
838 return changed_parameters
840 def check_state(self, atoms, tol=1e-15):
841 """Check for any system changes since last calculation."""
842 return compare_atoms(
843 self.atoms,
844 atoms,
845 tol=tol,
846 excluded_properties=set(self.ignored_changes),
847 )
849 def calculate(
850 self, atoms=None, properties=['energy'], system_changes=all_changes
851 ):
852 """Do the calculation.
854 properties: list of str
855 List of what needs to be calculated. Can be any combination
856 of 'energy', 'forces', 'stress', 'dipole', 'charges', 'magmom'
857 and 'magmoms'.
858 system_changes: list of str
859 List of what has changed since last calculation. Can be
860 any combination of these six: 'positions', 'numbers', 'cell',
861 'pbc', 'initial_charges' and 'initial_magmoms'.
863 Subclasses need to implement this, but can ignore properties
864 and system_changes if they want. Calculated properties should
865 be inserted into results dictionary like shown in this dummy
866 example::
868 self.results = {'energy': 0.0,
869 'forces': np.zeros((len(atoms), 3)),
870 'stress': np.zeros(6),
871 'dipole': np.zeros(3),
872 'charges': np.zeros(len(atoms)),
873 'magmom': 0.0,
874 'magmoms': np.zeros(len(atoms))}
876 The subclass implementation should first call this
877 implementation to set the atoms attribute and create any missing
878 directories.
879 """
880 if atoms is not None:
881 self.atoms = atoms.copy()
882 if not os.path.isdir(self._directory):
883 try:
884 os.makedirs(self._directory)
885 except FileExistsError as e:
886 # We can only end up here in case of a race condition if
887 # multiple Calculators are running concurrently *and* use the
888 # same _directory, which cannot be expected to work anyway.
889 msg = (
890 'Concurrent use of directory '
891 + self._directory
892 + 'by multiple Calculator instances detected. Please '
893 'use one directory per instance.'
894 )
895 raise RuntimeError(msg) from e
897 @deprecated('Please use `ase.calculators.fd.FiniteDifferenceCalculator`.')
898 def calculate_numerical_forces(self, atoms, d=0.001):
899 """Calculate numerical forces using finite difference.
901 All atoms will be displaced by +d and -d in all directions.
903 .. deprecated:: 3.24.0
904 """
905 from ase.calculators.fd import calculate_numerical_forces
907 return calculate_numerical_forces(atoms, eps=d)
909 @deprecated('Please use `ase.calculators.fd.FiniteDifferenceCalculator`.')
910 def calculate_numerical_stress(self, atoms, d=1e-6, voigt=True):
911 """Calculate numerical stress using finite difference.
913 .. deprecated:: 3.24.0
914 """
915 from ase.calculators.fd import calculate_numerical_stress
917 return calculate_numerical_stress(atoms, eps=d, voigt=voigt)
919 def _deprecated_get_spin_polarized(self):
920 msg = (
921 'This calculator does not implement get_spin_polarized(). '
922 'In the future, calc.get_spin_polarized() will work only on '
923 'calculator classes that explicitly implement this method or '
924 'inherit the method via specialized subclasses.'
925 )
926 warnings.warn(msg, FutureWarning)
927 return False
929 def band_structure(self):
930 """Create band-structure object for plotting."""
931 from ase.spectrum.band_structure import get_band_structure
933 # XXX This calculator is supposed to just have done a band structure
934 # calculation, but the calculator may not have the correct Fermi level
935 # if it updated the Fermi level after changing k-points.
936 # This will be a problem with some calculators (currently GPAW), and
937 # the user would have to override this by providing the Fermi level
938 # from the selfconsistent calculation.
939 return get_band_structure(calc=self)
942class OldShellProfile:
943 def __init__(self, command):
944 self.command = command
945 self.configvars = {}
947 def execute(self, calc):
948 if self.command is None:
949 raise EnvironmentError(
950 'Please set ${} environment variable '.format(
951 'ASE_' + self.calc.upper() + '_COMMAND'
952 )
953 + 'or supply the command keyword'
954 )
955 command = self.command
956 if 'PREFIX' in command:
957 command = command.replace('PREFIX', calc.prefix)
959 try:
960 proc = subprocess.Popen(command, shell=True, cwd=calc.directory)
961 except OSError as err:
962 # Actually this may never happen with shell=True, since
963 # probably the shell launches successfully. But we soon want
964 # to allow calling the subprocess directly, and then this
965 # distinction (failed to launch vs failed to run) is useful.
966 msg = f'Failed to execute "{command}"'
967 raise EnvironmentError(msg) from err
969 errorcode = proc.wait()
971 if errorcode:
972 path = os.path.abspath(calc.directory)
973 msg = (
974 'Calculator "{}" failed with command "{}" failed in '
975 '{} with error code {}'.format(
976 calc.name, command, path, errorcode
977 )
978 )
979 raise CalculationFailed(msg)
982@dataclass
983class FileIORules:
984 """Rules for controlling streams options to external command.
986 FileIOCalculator will direct stdin and stdout and append arguments
987 to the calculator command using the specifications on this class.
989 Currently names can contain "{prefix}" which will be substituted by
990 calc.prefix. This will go away if/when we can remove prefix."""
991 extend_argv: Sequence[str] = tuple()
992 stdin_name: Optional[str] = None
993 stdout_name: Optional[str] = None
995 configspec: Dict[str, Any] = field(default_factory=dict)
997 def load_config(self, section):
998 dct = {}
999 for key, value in self.configspec.items():
1000 if key in section:
1001 value = section[key]
1002 dct[key] = value
1003 return dct
1006class BadConfiguration(Exception):
1007 pass
1010def _validate_command(command: str) -> str:
1011 # We like to store commands as strings (and call shlex.split() later),
1012 # but we also like to validate them early. This will error out if
1013 # command contains syntax problems and will also normalize e.g.
1014 # multiple spaces:
1015 try:
1016 return shlex.join(shlex.split(command))
1017 except ValueError as err:
1018 raise BadConfiguration('Cannot parse command string') from err
1021@dataclass
1022class StandardProfile:
1023 command: str
1024 configvars: Dict[str, Any] = field(default_factory=dict)
1026 def __post_init__(self):
1027 self.command = _validate_command(self.command)
1029 def execute(self, calc):
1030 try:
1031 self._call(calc, subprocess.check_call)
1032 except subprocess.CalledProcessError as err:
1033 directory = Path(calc.directory).resolve()
1034 msg = (f'Calculator {calc.name} failed with args {err.args} '
1035 f'in directory {directory}')
1036 raise CalculationFailed(msg) from err
1038 def execute_nonblocking(self, calc):
1039 return self._call(calc, subprocess.Popen)
1041 @property
1042 def _split_command(self):
1043 # XXX Unduplicate common stuff between StandardProfile and
1044 # that of GenericFileIO
1045 return shlex.split(self.command)
1047 def _call(self, calc, subprocess_function):
1048 from contextlib import ExitStack
1050 directory = Path(calc.directory).resolve()
1051 fileio_rules = calc.fileio_rules
1053 with ExitStack() as stack:
1055 def _maybe_open(name, mode):
1056 if name is None:
1057 return None
1059 name = name.format(prefix=calc.prefix)
1060 directory = Path(calc.directory)
1061 return stack.enter_context(open(directory / name, mode))
1063 stdout_fd = _maybe_open(fileio_rules.stdout_name, 'wb')
1064 stdin_fd = _maybe_open(fileio_rules.stdin_name, 'rb')
1066 argv = [*self._split_command, *fileio_rules.extend_argv]
1067 argv = [arg.format(prefix=calc.prefix) for arg in argv]
1068 return subprocess_function(
1069 argv, cwd=directory,
1070 stdout=stdout_fd,
1071 stdin=stdin_fd)
1074class FileIOCalculator(Calculator):
1075 """Base class for calculators that write/read input/output files."""
1077 # Static specification of rules for this calculator:
1078 fileio_rules: Optional[FileIORules] = None
1080 # command: Optional[str] = None
1081 # 'Command used to start calculation'
1083 # Fallback command when nothing else is specified.
1084 # There will be no fallback in the future; it must be explicitly
1085 # configured.
1086 _legacy_default_command: Optional[str] = None
1088 cfg = _cfg # Ensure easy access to config for subclasses
1090 @classmethod
1091 def ruleset(cls, *args, **kwargs):
1092 """Helper for subclasses to define FileIORules."""
1093 return FileIORules(*args, **kwargs)
1095 def __init__(
1096 self,
1097 restart=None,
1098 ignore_bad_restart_file=Calculator._deprecated,
1099 label=None,
1100 atoms=None,
1101 command=None,
1102 profile=None,
1103 **kwargs,
1104 ):
1105 """File-IO calculator.
1107 command: str
1108 Command used to start calculation.
1109 """
1111 super().__init__(restart, ignore_bad_restart_file, label, atoms,
1112 **kwargs)
1114 if profile is None:
1115 profile = self._initialize_profile(command)
1116 self.profile = profile
1118 @property
1119 def command(self):
1120 # XXX deprecate me
1121 #
1122 # This is for calculators that invoke Popen directly on
1123 # self.command instead of letting us (superclass) do it.
1124 return self.profile.command
1126 @command.setter
1127 def command(self, command):
1128 self.profile.command = command
1130 @classmethod
1131 def load_argv_profile(cls, cfg, section_name):
1132 # Helper method to load configuration.
1133 # This is used by the tests, do not rely on this as it will change.
1134 try:
1135 section = cfg.parser[section_name]
1136 except KeyError:
1137 raise BadConfiguration(f'No {section_name!r} section')
1139 if cls.fileio_rules is not None:
1140 configvars = cls.fileio_rules.load_config(section)
1141 else:
1142 configvars = {}
1144 try:
1145 command = section['command']
1146 except KeyError:
1147 raise BadConfiguration(
1148 f'No command field in {section_name!r} section')
1150 return StandardProfile(command, configvars)
1152 def _initialize_profile(self, command):
1153 if command is None:
1154 name = 'ASE_' + self.name.upper() + '_COMMAND'
1155 command = self.cfg.get(name)
1157 if command is None and self.name in self.cfg.parser:
1158 return self.load_argv_profile(self.cfg, self.name)
1160 if command is None:
1161 # XXX issue a FutureWarning if this causes the command
1162 # to no longer be None
1163 command = self._legacy_default_command
1165 if command is None:
1166 raise EnvironmentError(
1167 f'No configuration of {self.name}. '
1168 f'Missing section [{self.name}] in configuration')
1170 return OldShellProfile(command)
1172 def calculate(
1173 self, atoms=None, properties=['energy'], system_changes=all_changes
1174 ):
1175 Calculator.calculate(self, atoms, properties, system_changes)
1176 self.write_input(self.atoms, properties, system_changes)
1177 self.execute()
1178 self.read_results()
1180 def execute(self):
1181 self.profile.execute(self)
1183 def write_input(self, atoms, properties=None, system_changes=None):
1184 """Write input file(s).
1186 Call this method first in subclasses so that directories are
1187 created automatically."""
1189 absdir = os.path.abspath(self.directory)
1190 if absdir != os.curdir and not os.path.isdir(self.directory):
1191 os.makedirs(self.directory)
1193 def read_results(self):
1194 """Read energy, forces, ... from output file(s)."""