Coverage for /builds/ase/ase/ase/calculators/calculator.py: 83.92%
541 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
1# fmt: off
3import copy
4import os
5import shlex
6import subprocess
7import warnings
8from abc import abstractmethod
9from dataclasses import dataclass, field
10from math import pi, sqrt
11from pathlib import Path
12from typing import Any, Dict, List, Optional, Sequence, Set, Union
14import numpy as np
16from ase.calculators.abc import GetPropertiesMixin
17from ase.cell import Cell
18from ase.config import cfg as _cfg
19from ase.outputs import Properties, all_outputs
20from ase.utils import deprecated, jsonable
22from .names import names
25class CalculatorError(RuntimeError):
26 """Base class of error types related to ASE calculators."""
29class CalculatorSetupError(CalculatorError):
30 """Calculation cannot be performed with the given parameters.
32 Reasons to raise this error are:
33 * The calculator is not properly configured
34 (missing executable, environment variables, ...)
35 * The given atoms object is not supported
36 * Calculator parameters are unsupported
38 Typically raised before a calculation."""
41class EnvironmentError(CalculatorSetupError):
42 """Raised if calculator is not properly set up with ASE.
44 May be missing an executable or environment variables."""
47class InputError(CalculatorSetupError):
48 """Raised if inputs given to the calculator were incorrect.
50 Bad input keywords or values, or missing pseudopotentials.
52 This may be raised before or during calculation, depending on
53 when the problem is detected."""
56class CalculationFailed(CalculatorError):
57 """Calculation failed unexpectedly.
59 Reasons to raise this error are:
60 * Calculation did not converge
61 * Calculation ran out of memory
62 * Segmentation fault or other abnormal termination
63 * Arithmetic trouble (singular matrices, NaN, ...)
65 Typically raised during calculation."""
68class SCFError(CalculationFailed):
69 """SCF loop did not converge."""
72class ReadError(CalculatorError):
73 """Unexpected irrecoverable error while reading calculation results."""
76class PropertyNotImplementedError(NotImplementedError):
77 """Raised if a calculator does not implement the requested property."""
80class PropertyNotPresent(CalculatorError):
81 """Requested property is missing.
83 Maybe it was never calculated, or for some reason was not extracted
84 with the rest of the results, without being a fatal ReadError."""
87def compare_atoms(atoms1, atoms2, tol=1e-15, excluded_properties=None):
88 """Check for system changes since last calculation. Properties in
89 ``excluded_properties`` are not checked."""
90 if atoms1 is None:
91 system_changes = all_changes[:]
92 else:
93 system_changes = []
95 properties_to_check = set(all_changes)
96 if excluded_properties:
97 properties_to_check -= set(excluded_properties)
99 # Check properties that aren't in Atoms.arrays but are attributes of
100 # Atoms objects
101 for prop in ['cell', 'pbc']:
102 if prop in properties_to_check:
103 properties_to_check.remove(prop)
104 if not equal(
105 getattr(atoms1, prop), getattr(atoms2, prop), atol=tol
106 ):
107 system_changes.append(prop)
109 arrays1 = set(atoms1.arrays)
110 arrays2 = set(atoms2.arrays)
112 # Add any properties that are only in atoms1.arrays or only in
113 # atoms2.arrays (and aren't excluded). Note that if, e.g. arrays1 has
114 # `initial_charges` which is merely zeros and arrays2 does not have
115 # this array, we'll still assume that the system has changed. However,
116 # this should only occur rarely.
117 system_changes += properties_to_check & (arrays1 ^ arrays2)
119 # Finally, check all of the non-excluded properties shared by the atoms
120 # arrays.
121 for prop in properties_to_check & arrays1 & arrays2:
122 if not equal(atoms1.arrays[prop], atoms2.arrays[prop], atol=tol):
123 system_changes.append(prop)
125 return system_changes
128all_properties = [
129 'energy',
130 'forces',
131 'stress',
132 'stresses',
133 'dipole',
134 'charges',
135 'magmom',
136 'magmoms',
137 'free_energy',
138 'energies',
139 'dielectric_tensor',
140 'born_effective_charges',
141 'polarization',
142]
145all_changes = [
146 'positions',
147 'numbers',
148 'cell',
149 'pbc',
150 'initial_charges',
151 'initial_magmoms',
152]
155special = {
156 'cp2k': 'CP2K',
157 'demonnano': 'DemonNano',
158 'dftd3': 'DFTD3',
159 'dmol': 'DMol3',
160 'eam': 'EAM',
161 'elk': 'ELK',
162 'emt': 'EMT',
163 'exciting': 'ExcitingGroundStateCalculator',
164 'crystal': 'CRYSTAL',
165 'ff': 'ForceField',
166 'gamess_us': 'GAMESSUS',
167 'gulp': 'GULP',
168 'kim': 'KIM',
169 'lammpsrun': 'LAMMPS',
170 'lammpslib': 'LAMMPSlib',
171 'lj': 'LennardJones',
172 'mopac': 'MOPAC',
173 'morse': 'MorsePotential',
174 'nwchem': 'NWChem',
175 'openmx': 'OpenMX',
176 'orca': 'ORCA',
177 'qchem': 'QChem',
178 'tip3p': 'TIP3P',
179 'tip4p': 'TIP4P',
180}
183external_calculators = {}
186def register_calculator_class(name, cls):
187 """Add the class into the database."""
188 assert name not in external_calculators
189 external_calculators[name] = cls
190 names.append(name)
191 names.sort()
194def get_calculator_class(name):
195 """Return calculator class."""
196 if name == 'asap':
197 from asap3 import EMT as Calculator
198 elif name == 'gpaw':
199 from gpaw import GPAW as Calculator
200 elif name == 'hotbit':
201 from hotbit import Calculator
202 elif name == 'vasp2':
203 from ase.calculators.vasp import Vasp2 as Calculator
204 elif name == 'ace':
205 from ase.calculators.acemolecule import ACE as Calculator
206 elif name == 'Psi4':
207 from ase.calculators.psi4 import Psi4 as Calculator
208 elif name == 'mattersim':
209 from mattersim.forcefield import MatterSimCalculator as Calculator
210 elif name == 'mace_mp':
211 from mace.calculators import mace_mp as Calculator
212 elif name in external_calculators:
213 Calculator = external_calculators[name]
214 else:
215 classname = special.get(name, name.title())
216 module = __import__('ase.calculators.' + name, {}, None, [classname])
217 Calculator = getattr(module, classname)
218 return Calculator
221def equal(a, b, tol=None, rtol=None, atol=None):
222 """ndarray-enabled comparison function."""
223 # XXX Known bugs:
224 # * Comparing cell objects (pbc not part of array representation)
225 # * Infinite recursion for cyclic dicts
226 # * Can of worms is open
227 if tol is not None:
228 msg = 'Use `equal(a, b, rtol=..., atol=...)` instead of `tol=...`'
229 warnings.warn(msg, DeprecationWarning)
230 assert (
231 rtol is None and atol is None
232 ), 'Do not use deprecated `tol` with `atol` and/or `rtol`'
233 rtol = tol
234 atol = tol
236 a_is_dict = isinstance(a, dict)
237 b_is_dict = isinstance(b, dict)
238 if a_is_dict or b_is_dict:
239 # Check that both a and b are dicts
240 if not (a_is_dict and b_is_dict):
241 return False
242 if a.keys() != b.keys():
243 return False
244 return all(equal(a[key], b[key], rtol=rtol, atol=atol) for key in a)
246 if np.shape(a) != np.shape(b):
247 return False
249 if rtol is None and atol is None:
250 return np.array_equal(a, b)
252 if rtol is None:
253 rtol = 0
254 if atol is None:
255 atol = 0
257 return np.allclose(a, b, rtol=rtol, atol=atol)
260def kptdensity2monkhorstpack(atoms, kptdensity=3.5, even=True):
261 """Convert k-point density to Monkhorst-Pack grid size.
263 atoms: Atoms object
264 Contains unit cell and information about boundary conditions.
265 kptdensity: float
266 Required k-point density. Default value is 3.5 point per Ang^-1.
267 even: bool
268 Round up to even numbers.
269 """
271 recipcell = atoms.cell.reciprocal()
272 kpts = []
273 for i in range(3):
274 if atoms.pbc[i]:
275 k = 2 * pi * sqrt((recipcell[i] ** 2).sum()) * kptdensity
276 if even:
277 kpts.append(2 * int(np.ceil(k / 2)))
278 else:
279 kpts.append(int(np.ceil(k)))
280 else:
281 kpts.append(1)
282 return np.array(kpts)
285def kpts2mp(atoms, kpts, even=False):
286 if kpts is None:
287 return np.array([1, 1, 1])
288 if isinstance(kpts, (float, int)):
289 return kptdensity2monkhorstpack(atoms, kpts, even)
290 else:
291 return kpts
294def kpts2sizeandoffsets(
295 size=None, density=None, gamma=None, even=None, atoms=None
296):
297 """Helper function for selecting k-points.
299 Use either size or density.
301 size: 3 ints
302 Number of k-points.
303 density: float
304 K-point density in units of k-points per Ang^-1.
305 gamma: None or bool
306 Should the Gamma-point be included? Yes / no / don't care:
307 True / False / None.
308 even: None or bool
309 Should the number of k-points be even? Yes / no / don't care:
310 True / False / None.
311 atoms: Atoms object
312 Needed for calculating k-point density.
314 """
316 if size is not None and density is not None:
317 raise ValueError(
318 'Cannot specify k-point mesh size and ' 'density simultaneously'
319 )
320 elif density is not None and atoms is None:
321 raise ValueError(
322 'Cannot set k-points from "density" unless '
323 'Atoms are provided (need BZ dimensions).'
324 )
326 if size is None:
327 if density is None:
328 size = [1, 1, 1]
329 else:
330 size = kptdensity2monkhorstpack(atoms, density, None)
332 # Not using the rounding from kptdensity2monkhorstpack as it doesn't do
333 # rounding to odd numbers
334 if even is not None:
335 size = np.array(size)
336 remainder = size % 2
337 if even:
338 size += remainder
339 else: # Round up to odd numbers
340 size += 1 - remainder
342 offsets = [0, 0, 0]
343 if atoms is None:
344 pbc = [True, True, True]
345 else:
346 pbc = atoms.pbc
348 if gamma is not None:
349 for i, s in enumerate(size):
350 if pbc[i] and s % 2 != bool(gamma):
351 offsets[i] = 0.5 / s
353 return size, offsets
356@jsonable('kpoints')
357class KPoints:
358 def __init__(self, kpts=None):
359 if kpts is None:
360 kpts = np.zeros((1, 3))
361 self.kpts = kpts
363 def todict(self):
364 return vars(self)
367def kpts2kpts(kpts, atoms=None):
368 from ase.dft.kpoints import monkhorst_pack
370 if kpts is None:
371 return KPoints()
373 if hasattr(kpts, 'kpts'):
374 return kpts
376 if isinstance(kpts, dict):
377 if 'kpts' in kpts:
378 return KPoints(kpts['kpts'])
379 if 'path' in kpts:
380 cell = Cell.ascell(atoms.cell)
381 return cell.bandpath(pbc=atoms.pbc, **kpts)
382 size, offsets = kpts2sizeandoffsets(atoms=atoms, **kpts)
383 return KPoints(monkhorst_pack(size) + offsets)
385 if isinstance(kpts[0], int):
386 return KPoints(monkhorst_pack(kpts))
388 return KPoints(np.array(kpts))
391def kpts2ndarray(kpts, atoms=None):
392 """Convert kpts keyword to 2-d ndarray of scaled k-points."""
393 return kpts2kpts(kpts, atoms=atoms).kpts
396class Parameters(dict):
397 """Dictionary for parameters.
399 Special feature: If param is a Parameters instance, then param.xc
400 is a shorthand for param['xc'].
401 """
403 def __getattr__(self, key):
404 if key not in self:
405 return dict.__getattribute__(self, key)
406 return self[key]
408 def __setattr__(self, key, value):
409 self[key] = value
411 @classmethod
412 def read(cls, filename):
413 """Read parameters from file."""
414 # We use ast to evaluate literals, avoiding eval()
415 # for security reasons.
416 import ast
418 with open(filename) as fd:
419 txt = fd.read().strip()
420 assert txt.startswith('dict(')
421 assert txt.endswith(')')
422 txt = txt[5:-1]
424 # The tostring() representation "dict(...)" is not actually
425 # a literal, so we manually parse that along with the other
426 # formatting that we did manually:
427 dct = {}
428 for line in txt.splitlines():
429 key, val = line.split('=', 1)
430 key = key.strip()
431 val = val.strip()
432 if val[-1] == ',':
433 val = val[:-1]
434 dct[key] = ast.literal_eval(val)
436 parameters = cls(dct)
437 return parameters
439 def tostring(self):
440 keys = sorted(self)
441 return (
442 'dict('
443 + ',\n '.join(f'{key}={self[key]!r}' for key in keys)
444 + ')\n'
445 )
447 def write(self, filename):
448 Path(filename).write_text(self.tostring())
451class BaseCalculator(GetPropertiesMixin):
452 implemented_properties: List[str] = []
453 'Properties calculator can handle (energy, forces, ...)'
455 # Placeholder object for deprecated arguments. Let deprecated keywords
456 # default to _deprecated and then issue a warning if the user passed
457 # any other object (such as None).
458 _deprecated = object()
460 def __init__(self, parameters=None, use_cache=True):
461 if parameters is None:
462 parameters = {}
464 self.parameters = dict(parameters)
465 self.atoms = None
466 self.results = {}
467 self.use_cache = use_cache
469 def calculate_properties(self, atoms, properties):
470 """This method is experimental; currently for internal use."""
471 for name in properties:
472 if name not in all_outputs:
473 raise ValueError(f'No such property: {name}')
475 # We ignore system changes for now.
476 self.calculate(atoms, properties, system_changes=all_changes)
478 props = self.export_properties()
480 for name in properties:
481 if name not in props:
482 raise PropertyNotPresent(name)
483 return props
485 @abstractmethod
486 def calculate(self, atoms, properties, system_changes):
487 ...
489 def check_state(self, atoms, tol=1e-15):
490 """Check for any system changes since last calculation."""
491 if self.use_cache:
492 return compare_atoms(self.atoms, atoms, tol=tol)
493 else:
494 return all_changes
496 def get_property(self, name, atoms=None, allow_calculation=True):
497 if name not in self.implemented_properties:
498 raise PropertyNotImplementedError(
499 f'{name} property not implemented'
500 )
502 if atoms is None:
503 atoms = self.atoms
504 system_changes = []
505 else:
506 system_changes = self.check_state(atoms)
508 if system_changes:
509 self.atoms = None
510 self.results = {}
512 if name not in self.results:
513 if not allow_calculation:
514 return None
516 if self.use_cache:
517 self.atoms = atoms.copy()
519 self.calculate(atoms, [name], system_changes)
521 if name not in self.results:
522 # For some reason the calculator was not able to do what we want,
523 # and that is OK.
524 raise PropertyNotImplementedError(
525 '{} not present in this ' 'calculation'.format(name)
526 )
528 result = self.results[name]
529 if isinstance(result, np.ndarray):
530 result = result.copy()
531 return result
533 def calculation_required(self, atoms, properties):
534 assert not isinstance(properties, str)
535 system_changes = self.check_state(atoms)
536 if system_changes:
537 return True
538 for name in properties:
539 if name not in self.results:
540 return True
541 return False
543 def export_properties(self):
544 return Properties(self.results)
546 def _get_name(self) -> str: # child class can override this
547 return self.__class__.__name__.lower()
549 @property
550 def name(self) -> str:
551 return self._get_name()
553 def todict(self) -> dict[str, Any]:
554 """Obtain a dictionary of parameter information"""
555 return {}
558class Calculator(BaseCalculator):
559 """Base-class for all ASE calculators.
561 A calculator must raise PropertyNotImplementedError if asked for a
562 property that it can't calculate. So, if calculation of the
563 stress tensor has not been implemented, get_stress(atoms) should
564 raise PropertyNotImplementedError. This can be achieved simply by not
565 including the string 'stress' in the list implemented_properties
566 which is a class member. These are the names of the standard
567 properties: 'energy', 'forces', 'stress', 'dipole', 'charges',
568 'magmom' and 'magmoms'.
569 """
571 default_parameters: Dict[str, Any] = {}
572 'Default parameters'
574 ignored_changes: Set[str] = set()
575 'Properties of Atoms which we ignore for the purposes of cache '
576 'invalidation with check_state().'
578 discard_results_on_any_change = False
579 'Whether we purge the results following any change in the set() method. '
580 'Most (file I/O) calculators will probably want this.'
582 def __init__(
583 self,
584 restart=None,
585 ignore_bad_restart_file=BaseCalculator._deprecated,
586 label=None,
587 atoms=None,
588 directory='.',
589 **kwargs,
590 ):
591 """Basic calculator implementation.
593 restart: str
594 Prefix for restart file. May contain a directory. Default
595 is None: don't restart.
596 ignore_bad_restart_file: bool
597 Deprecated, please do not use.
598 Passing more than one positional argument to Calculator()
599 is deprecated and will stop working in the future.
600 Ignore broken or missing restart file. By default, it is an
601 error if the restart file is missing or broken.
602 directory: str or PurePath
603 Working directory in which to read and write files and
604 perform calculations.
605 label: str
606 Name used for all files. Not supported by all calculators.
607 May contain a directory, but please use the directory parameter
608 for that instead.
609 atoms: Atoms object
610 Optional Atoms object to which the calculator will be
611 attached. When restarting, atoms will get its positions and
612 unit-cell updated from file.
613 """
614 self.atoms = None # copy of atoms object from last calculation
615 self.results = {} # calculated properties (energy, forces, ...)
616 self.parameters = None # calculational parameters
617 self._directory = None # Initialize
619 if ignore_bad_restart_file is self._deprecated:
620 ignore_bad_restart_file = False
621 else:
622 warnings.warn(
623 FutureWarning(
624 'The keyword "ignore_bad_restart_file" is deprecated and '
625 'will be removed in a future version of ASE. Passing more '
626 'than one positional argument to Calculator is also '
627 'deprecated and will stop functioning in the future. '
628 'Please pass arguments by keyword (key=value) except '
629 'optionally the "restart" keyword.'
630 )
631 )
633 if restart is not None:
634 try:
635 self.read(restart) # read parameters, atoms and results
636 except ReadError:
637 if ignore_bad_restart_file:
638 self.reset()
639 else:
640 raise
642 self.directory = directory
643 self.prefix = None
644 if label is not None:
645 if self.directory == '.' and '/' in label:
646 # We specified directory in label, and nothing in the diretory
647 # key
648 self.label = label
649 elif '/' not in label:
650 # We specified our directory in the directory keyword
651 # or not at all
652 self.label = '/'.join((self.directory, label))
653 else:
654 raise ValueError(
655 'Directory redundantly specified though '
656 'directory="{}" and label="{}". '
657 'Please omit "/" in label.'.format(self.directory, label)
658 )
660 if self.parameters is None:
661 # Use default parameters if they were not read from file:
662 self.parameters = self.get_default_parameters()
664 if atoms is not None:
665 atoms.calc = self
666 if self.atoms is not None:
667 # Atoms were read from file. Update atoms:
668 if not (
669 equal(atoms.numbers, self.atoms.numbers)
670 and (atoms.pbc == self.atoms.pbc).all()
671 ):
672 raise CalculatorError('Atoms not compatible with file')
673 atoms.positions = self.atoms.positions
674 atoms.cell = self.atoms.cell
676 self.set(**kwargs)
678 if not hasattr(self, 'get_spin_polarized'):
679 self.get_spin_polarized = self._deprecated_get_spin_polarized
680 # XXX We are very naughty and do not call super constructor!
682 # For historical reasons we have a particular caching protocol.
683 # We disable the superclass' optional cache.
684 self.use_cache = False
686 @property
687 def directory(self) -> str:
688 return self._directory
690 @directory.setter
691 def directory(self, directory: Union[str, os.PathLike]):
692 self._directory = str(Path(directory)) # Normalize path.
694 @property
695 def label(self):
696 if self.directory == '.':
697 return self.prefix
699 # Generally, label ~ directory/prefix
700 #
701 # We use '/' rather than os.pathsep because
702 # 1) directory/prefix does not represent any actual path
703 # 2) We want the same string to work the same on all platforms
704 if self.prefix is None:
705 return self.directory + '/'
707 return f'{self.directory}/{self.prefix}'
709 @label.setter
710 def label(self, label):
711 if label is None:
712 self.directory = '.'
713 self.prefix = None
714 return
716 tokens = label.rsplit('/', 1)
717 if len(tokens) == 2:
718 directory, prefix = tokens
719 else:
720 assert len(tokens) == 1
721 directory = '.'
722 prefix = tokens[0]
723 if prefix == '':
724 prefix = None
725 self.directory = directory
726 self.prefix = prefix
728 def set_label(self, label):
729 """Set label and convert label to directory and prefix.
731 Examples:
733 * label='abc': (directory='.', prefix='abc')
734 * label='dir1/abc': (directory='dir1', prefix='abc')
735 * label=None: (directory='.', prefix=None)
736 """
737 self.label = label
739 def get_default_parameters(self):
740 return Parameters(copy.deepcopy(self.default_parameters))
742 def todict(self, skip_default=True):
743 defaults = self.get_default_parameters()
744 dct = {}
745 for key, value in self.parameters.items():
746 if hasattr(value, 'todict'):
747 value = value.todict()
748 if skip_default:
749 default = defaults.get(key, '_no_default_')
750 if default != '_no_default_' and equal(value, default):
751 continue
752 dct[key] = value
753 return dct
755 def reset(self):
756 """Clear all information from old calculation."""
758 self.atoms = None
759 self.results = {}
761 def read(self, label):
762 """Read atoms, parameters and calculated properties from output file.
764 Read result from self.label file. Raise ReadError if the file
765 is not there. If the file is corrupted or contains an error
766 message from the calculation, a ReadError should also be
767 raised. In case of succes, these attributes must set:
769 atoms: Atoms object
770 The state of the atoms from last calculation.
771 parameters: Parameters object
772 The parameter dictionary.
773 results: dict
774 Calculated properties like energy and forces.
776 The FileIOCalculator.read() method will typically read atoms
777 and parameters and get the results dict by calling the
778 read_results() method."""
780 self.set_label(label)
782 def get_atoms(self):
783 if self.atoms is None:
784 raise ValueError('Calculator has no atoms')
785 atoms = self.atoms.copy()
786 atoms.calc = self
787 return atoms
789 @classmethod
790 def read_atoms(cls, restart, **kwargs):
791 return cls(restart=restart, label=restart, **kwargs).get_atoms()
793 def set(self, **kwargs):
794 """Set parameters like set(key1=value1, key2=value2, ...).
796 A dictionary containing the parameters that have been changed
797 is returned.
799 Subclasses must implement a set() method that will look at the
800 chaneged parameters and decide if a call to reset() is needed.
801 If the changed parameters are harmless, like a change in
802 verbosity, then there is no need to call reset().
804 The special keyword 'parameters' can be used to read
805 parameters from a file."""
807 if 'parameters' in kwargs:
808 filename = kwargs.pop('parameters')
809 parameters = Parameters.read(filename)
810 parameters.update(kwargs)
811 kwargs = parameters
813 changed_parameters = {}
815 for key, value in kwargs.items():
816 oldvalue = self.parameters.get(key)
817 if key not in self.parameters or not equal(value, oldvalue):
818 changed_parameters[key] = value
819 self.parameters[key] = value
821 if self.discard_results_on_any_change and changed_parameters:
822 self.reset()
823 return changed_parameters
825 def check_state(self, atoms, tol=1e-15):
826 """Check for any system changes since last calculation."""
827 return compare_atoms(
828 self.atoms,
829 atoms,
830 tol=tol,
831 excluded_properties=set(self.ignored_changes),
832 )
834 def calculate(
835 self, atoms=None, properties=['energy'], system_changes=all_changes
836 ):
837 """Do the calculation.
839 properties: list of str
840 List of what needs to be calculated. Can be any combination
841 of 'energy', 'forces', 'stress', 'dipole', 'charges', 'magmom'
842 and 'magmoms'.
843 system_changes: list of str
844 List of what has changed since last calculation. Can be
845 any combination of these six: 'positions', 'numbers', 'cell',
846 'pbc', 'initial_charges' and 'initial_magmoms'.
848 Subclasses need to implement this, but can ignore properties
849 and system_changes if they want. Calculated properties should
850 be inserted into results dictionary like shown in this dummy
851 example::
853 self.results = {'energy': 0.0,
854 'forces': np.zeros((len(atoms), 3)),
855 'stress': np.zeros(6),
856 'dipole': np.zeros(3),
857 'charges': np.zeros(len(atoms)),
858 'magmom': 0.0,
859 'magmoms': np.zeros(len(atoms))}
861 The subclass implementation should first call this
862 implementation to set the atoms attribute and create any missing
863 directories.
864 """
865 if atoms is not None:
866 self.atoms = atoms.copy()
867 if not os.path.isdir(self._directory):
868 try:
869 os.makedirs(self._directory)
870 except FileExistsError as e:
871 # We can only end up here in case of a race condition if
872 # multiple Calculators are running concurrently *and* use the
873 # same _directory, which cannot be expected to work anyway.
874 msg = (
875 'Concurrent use of directory '
876 + self._directory
877 + 'by multiple Calculator instances detected. Please '
878 'use one directory per instance.'
879 )
880 raise RuntimeError(msg) from e
882 @deprecated('Please use `ase.calculators.fd.FiniteDifferenceCalculator`.')
883 def calculate_numerical_forces(self, atoms, d=0.001):
884 """Calculate numerical forces using finite difference.
886 All atoms will be displaced by +d and -d in all directions.
888 .. deprecated:: 3.24.0
889 """
890 from ase.calculators.fd import calculate_numerical_forces
892 return calculate_numerical_forces(atoms, eps=d)
894 @deprecated('Please use `ase.calculators.fd.FiniteDifferenceCalculator`.')
895 def calculate_numerical_stress(self, atoms, d=1e-6, voigt=True):
896 """Calculate numerical stress using finite difference.
898 .. deprecated:: 3.24.0
899 """
900 from ase.calculators.fd import calculate_numerical_stress
902 return calculate_numerical_stress(atoms, eps=d, voigt=voigt)
904 def _deprecated_get_spin_polarized(self):
905 msg = (
906 'This calculator does not implement get_spin_polarized(). '
907 'In the future, calc.get_spin_polarized() will work only on '
908 'calculator classes that explicitly implement this method or '
909 'inherit the method via specialized subclasses.'
910 )
911 warnings.warn(msg, FutureWarning)
912 return False
914 def band_structure(self):
915 """Create band-structure object for plotting."""
916 from ase.spectrum.band_structure import get_band_structure
918 # XXX This calculator is supposed to just have done a band structure
919 # calculation, but the calculator may not have the correct Fermi level
920 # if it updated the Fermi level after changing k-points.
921 # This will be a problem with some calculators (currently GPAW), and
922 # the user would have to override this by providing the Fermi level
923 # from the selfconsistent calculation.
924 return get_band_structure(calc=self)
927class OldShellProfile:
928 def __init__(self, command):
929 self.command = command
930 self.configvars = {}
932 def execute(self, calc):
933 if self.command is None:
934 raise EnvironmentError(
935 'Please set ${} environment variable '.format(
936 'ASE_' + self.calc.upper() + '_COMMAND'
937 )
938 + 'or supply the command keyword'
939 )
940 command = self.command
941 if 'PREFIX' in command:
942 command = command.replace('PREFIX', calc.prefix)
944 try:
945 proc = subprocess.Popen(command, shell=True, cwd=calc.directory)
946 except OSError as err:
947 # Actually this may never happen with shell=True, since
948 # probably the shell launches successfully. But we soon want
949 # to allow calling the subprocess directly, and then this
950 # distinction (failed to launch vs failed to run) is useful.
951 msg = f'Failed to execute "{command}"'
952 raise EnvironmentError(msg) from err
954 errorcode = proc.wait()
956 if errorcode:
957 path = os.path.abspath(calc.directory)
958 msg = (
959 'Calculator "{}" failed with command "{}" failed in '
960 '{} with error code {}'.format(
961 calc.name, command, path, errorcode
962 )
963 )
964 raise CalculationFailed(msg)
967@dataclass
968class FileIORules:
969 """Rules for controlling streams options to external command.
971 FileIOCalculator will direct stdin and stdout and append arguments
972 to the calculator command using the specifications on this class.
974 Currently names can contain "{prefix}" which will be substituted by
975 calc.prefix. This will go away if/when we can remove prefix."""
976 extend_argv: Sequence[str] = tuple()
977 stdin_name: Optional[str] = None
978 stdout_name: Optional[str] = None
980 configspec: Dict[str, Any] = field(default_factory=dict)
982 def load_config(self, section):
983 dct = {}
984 for key, value in self.configspec.items():
985 if key in section:
986 value = section[key]
987 dct[key] = value
988 return dct
991class BadConfiguration(Exception):
992 pass
995def _validate_command(command: str) -> str:
996 # We like to store commands as strings (and call shlex.split() later),
997 # but we also like to validate them early. This will error out if
998 # command contains syntax problems and will also normalize e.g.
999 # multiple spaces:
1000 try:
1001 return shlex.join(shlex.split(command))
1002 except ValueError as err:
1003 raise BadConfiguration('Cannot parse command string') from err
1006@dataclass
1007class StandardProfile:
1008 command: str
1009 configvars: Dict[str, Any] = field(default_factory=dict)
1011 def __post_init__(self):
1012 self.command = _validate_command(self.command)
1014 def execute(self, calc):
1015 try:
1016 self._call(calc, subprocess.check_call)
1017 except subprocess.CalledProcessError as err:
1018 directory = Path(calc.directory).resolve()
1019 msg = (f'Calculator {calc.name} failed with args {err.args} '
1020 f'in directory {directory}')
1021 raise CalculationFailed(msg) from err
1023 def execute_nonblocking(self, calc):
1024 return self._call(calc, subprocess.Popen)
1026 @property
1027 def _split_command(self):
1028 # XXX Unduplicate common stuff between StandardProfile and
1029 # that of GenericFileIO
1030 return shlex.split(self.command)
1032 def _call(self, calc, subprocess_function):
1033 from contextlib import ExitStack
1035 directory = Path(calc.directory).resolve()
1036 fileio_rules = calc.fileio_rules
1038 with ExitStack() as stack:
1040 def _maybe_open(name, mode):
1041 if name is None:
1042 return None
1044 name = name.format(prefix=calc.prefix)
1045 directory = Path(calc.directory)
1046 return stack.enter_context(open(directory / name, mode))
1048 stdout_fd = _maybe_open(fileio_rules.stdout_name, 'wb')
1049 stdin_fd = _maybe_open(fileio_rules.stdin_name, 'rb')
1051 argv = [*self._split_command, *fileio_rules.extend_argv]
1052 argv = [arg.format(prefix=calc.prefix) for arg in argv]
1053 return subprocess_function(
1054 argv, cwd=directory,
1055 stdout=stdout_fd,
1056 stdin=stdin_fd)
1059class FileIOCalculator(Calculator):
1060 """Base class for calculators that write/read input/output files."""
1062 # Static specification of rules for this calculator:
1063 fileio_rules: Optional[FileIORules] = None
1065 # command: Optional[str] = None
1066 # 'Command used to start calculation'
1068 # Fallback command when nothing else is specified.
1069 # There will be no fallback in the future; it must be explicitly
1070 # configured.
1071 _legacy_default_command: Optional[str] = None
1073 cfg = _cfg # Ensure easy access to config for subclasses
1075 @classmethod
1076 def ruleset(cls, *args, **kwargs):
1077 """Helper for subclasses to define FileIORules."""
1078 return FileIORules(*args, **kwargs)
1080 def __init__(
1081 self,
1082 restart=None,
1083 ignore_bad_restart_file=Calculator._deprecated,
1084 label=None,
1085 atoms=None,
1086 command=None,
1087 profile=None,
1088 **kwargs,
1089 ):
1090 """File-IO calculator.
1092 command: str
1093 Command used to start calculation.
1094 """
1096 super().__init__(restart, ignore_bad_restart_file, label, atoms,
1097 **kwargs)
1099 if profile is None:
1100 profile = self._initialize_profile(command)
1101 self.profile = profile
1103 @property
1104 def command(self):
1105 # XXX deprecate me
1106 #
1107 # This is for calculators that invoke Popen directly on
1108 # self.command instead of letting us (superclass) do it.
1109 return self.profile.command
1111 @command.setter
1112 def command(self, command):
1113 self.profile.command = command
1115 @classmethod
1116 def load_argv_profile(cls, cfg, section_name):
1117 # Helper method to load configuration.
1118 # This is used by the tests, do not rely on this as it will change.
1119 try:
1120 section = cfg.parser[section_name]
1121 except KeyError:
1122 raise BadConfiguration(f'No {section_name!r} section')
1124 if cls.fileio_rules is not None:
1125 configvars = cls.fileio_rules.load_config(section)
1126 else:
1127 configvars = {}
1129 try:
1130 command = section['command']
1131 except KeyError:
1132 raise BadConfiguration(
1133 f'No command field in {section_name!r} section')
1135 return StandardProfile(command, configvars)
1137 def _initialize_profile(self, command):
1138 if command is None:
1139 name = 'ASE_' + self.name.upper() + '_COMMAND'
1140 command = self.cfg.get(name)
1142 if command is None and self.name in self.cfg.parser:
1143 return self.load_argv_profile(self.cfg, self.name)
1145 if command is None:
1146 # XXX issue a FutureWarning if this causes the command
1147 # to no longer be None
1148 command = self._legacy_default_command
1150 if command is None:
1151 raise EnvironmentError(
1152 f'No configuration of {self.name}. '
1153 f'Missing section [{self.name}] in configuration')
1155 return OldShellProfile(command)
1157 def calculate(
1158 self, atoms=None, properties=['energy'], system_changes=all_changes
1159 ):
1160 Calculator.calculate(self, atoms, properties, system_changes)
1161 self.write_input(self.atoms, properties, system_changes)
1162 self.execute()
1163 self.read_results()
1165 def execute(self):
1166 self.profile.execute(self)
1168 def write_input(self, atoms, properties=None, system_changes=None):
1169 """Write input file(s).
1171 Call this method first in subclasses so that directories are
1172 created automatically."""
1174 absdir = os.path.abspath(self.directory)
1175 if absdir != os.curdir and not os.path.isdir(self.directory):
1176 os.makedirs(self.directory)
1178 def read_results(self):
1179 """Read energy, forces, ... from output file(s)."""