Coverage for /builds/kinetik161/ase/ase/calculators/calculator.py: 82.75%
487 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
1import copy
2import os
3import subprocess
4import warnings
6from abc import abstractmethod
7from math import pi, sqrt
8from pathlib import Path
9from typing import Any, Dict, List, Optional, Set, Union
11import numpy as np
13from ase.calculators.abc import GetPropertiesMixin
14from ase.cell import Cell
15from ase.config import cfg as _cfg
16from ase.outputs import Properties, all_outputs
17from ase.utils import jsonable
19from .names import names
22class CalculatorError(RuntimeError):
23 """Base class of error types related to ASE calculators."""
26class CalculatorSetupError(CalculatorError):
27 """Calculation cannot be performed with the given parameters.
29 Reasons to raise this error are:
30 * The calculator is not properly configured
31 (missing executable, environment variables, ...)
32 * The given atoms object is not supported
33 * Calculator parameters are unsupported
35 Typically raised before a calculation."""
38class EnvironmentError(CalculatorSetupError):
39 """Raised if calculator is not properly set up with ASE.
41 May be missing an executable or environment variables."""
44class InputError(CalculatorSetupError):
45 """Raised if inputs given to the calculator were incorrect.
47 Bad input keywords or values, or missing pseudopotentials.
49 This may be raised before or during calculation, depending on
50 when the problem is detected."""
53class CalculationFailed(CalculatorError):
54 """Calculation failed unexpectedly.
56 Reasons to raise this error are:
57 * Calculation did not converge
58 * Calculation ran out of memory
59 * Segmentation fault or other abnormal termination
60 * Arithmetic trouble (singular matrices, NaN, ...)
62 Typically raised during calculation."""
65class SCFError(CalculationFailed):
66 """SCF loop did not converge."""
69class ReadError(CalculatorError):
70 """Unexpected irrecoverable error while reading calculation results."""
73class PropertyNotImplementedError(NotImplementedError):
74 """Raised if a calculator does not implement the requested property."""
77class PropertyNotPresent(CalculatorError):
78 """Requested property is missing.
80 Maybe it was never calculated, or for some reason was not extracted
81 with the rest of the results, without being a fatal ReadError."""
84def compare_atoms(atoms1, atoms2, tol=1e-15, excluded_properties=None):
85 """Check for system changes since last calculation. Properties in
86 ``excluded_properties`` are not checked."""
87 if atoms1 is None:
88 system_changes = all_changes[:]
89 else:
90 system_changes = []
92 properties_to_check = set(all_changes)
93 if excluded_properties:
94 properties_to_check -= set(excluded_properties)
96 # Check properties that aren't in Atoms.arrays but are attributes of
97 # Atoms objects
98 for prop in ['cell', 'pbc']:
99 if prop in properties_to_check:
100 properties_to_check.remove(prop)
101 if not equal(
102 getattr(atoms1, prop), getattr(atoms2, prop), atol=tol
103 ):
104 system_changes.append(prop)
106 arrays1 = set(atoms1.arrays)
107 arrays2 = set(atoms2.arrays)
109 # Add any properties that are only in atoms1.arrays or only in
110 # atoms2.arrays (and aren't excluded). Note that if, e.g. arrays1 has
111 # `initial_charges` which is merely zeros and arrays2 does not have
112 # this array, we'll still assume that the system has changed. However,
113 # this should only occur rarely.
114 system_changes += properties_to_check & (arrays1 ^ arrays2)
116 # Finally, check all of the non-excluded properties shared by the atoms
117 # arrays.
118 for prop in properties_to_check & arrays1 & arrays2:
119 if not equal(atoms1.arrays[prop], atoms2.arrays[prop], atol=tol):
120 system_changes.append(prop)
122 return system_changes
125all_properties = [
126 'energy',
127 'forces',
128 'stress',
129 'stresses',
130 'dipole',
131 'charges',
132 'magmom',
133 'magmoms',
134 'free_energy',
135 'energies',
136 'dielectric_tensor',
137 'born_effective_charges',
138 'polarization',
139]
142all_changes = [
143 'positions',
144 'numbers',
145 'cell',
146 'pbc',
147 'initial_charges',
148 'initial_magmoms',
149]
152special = {
153 'cp2k': 'CP2K',
154 'demonnano': 'DemonNano',
155 'dftd3': 'DFTD3',
156 'dmol': 'DMol3',
157 'eam': 'EAM',
158 'elk': 'ELK',
159 'emt': 'EMT',
160 'exciting': 'ExcitingGroundStateCalculator',
161 'crystal': 'CRYSTAL',
162 'ff': 'ForceField',
163 'gamess_us': 'GAMESSUS',
164 'gulp': 'GULP',
165 'kim': 'KIM',
166 'lammpsrun': 'LAMMPS',
167 'lammpslib': 'LAMMPSlib',
168 'lj': 'LennardJones',
169 'mopac': 'MOPAC',
170 'morse': 'MorsePotential',
171 'nwchem': 'NWChem',
172 'openmx': 'OpenMX',
173 'orca': 'ORCA',
174 'qchem': 'QChem',
175 'tip3p': 'TIP3P',
176 'tip4p': 'TIP4P',
177}
180external_calculators = {}
183def register_calculator_class(name, cls):
184 """Add the class into the database."""
185 assert name not in external_calculators
186 external_calculators[name] = cls
187 names.append(name)
188 names.sort()
191def get_calculator_class(name):
192 """Return calculator class."""
193 if name == 'asap':
194 from asap3 import EMT as Calculator
195 elif name == 'gpaw':
196 from gpaw import GPAW as Calculator
197 elif name == 'hotbit':
198 from hotbit import Calculator
199 elif name == 'vasp2':
200 from ase.calculators.vasp import Vasp2 as Calculator
201 elif name == 'ace':
202 from ase.calculators.acemolecule import ACE as Calculator
203 elif name == 'Psi4':
204 from ase.calculators.psi4 import Psi4 as Calculator
205 elif name in external_calculators:
206 Calculator = external_calculators[name]
207 else:
208 classname = special.get(name, name.title())
209 module = __import__('ase.calculators.' + name, {}, None, [classname])
210 Calculator = getattr(module, classname)
211 return Calculator
214def equal(a, b, tol=None, rtol=None, atol=None):
215 """ndarray-enabled comparison function."""
216 # XXX Known bugs:
217 # * Comparing cell objects (pbc not part of array representation)
218 # * Infinite recursion for cyclic dicts
219 # * Can of worms is open
220 if tol is not None:
221 msg = 'Use `equal(a, b, rtol=..., atol=...)` instead of `tol=...`'
222 warnings.warn(msg, DeprecationWarning)
223 assert (
224 rtol is None and atol is None
225 ), 'Do not use deprecated `tol` with `atol` and/or `rtol`'
226 rtol = tol
227 atol = tol
229 a_is_dict = isinstance(a, dict)
230 b_is_dict = isinstance(b, dict)
231 if a_is_dict or b_is_dict:
232 # Check that both a and b are dicts
233 if not (a_is_dict and b_is_dict):
234 return False
235 if a.keys() != b.keys():
236 return False
237 return all(equal(a[key], b[key], rtol=rtol, atol=atol) for key in a)
239 if np.shape(a) != np.shape(b):
240 return False
242 if rtol is None and atol is None:
243 return np.array_equal(a, b)
245 if rtol is None:
246 rtol = 0
247 if atol is None:
248 atol = 0
250 return np.allclose(a, b, rtol=rtol, atol=atol)
253def kptdensity2monkhorstpack(atoms, kptdensity=3.5, even=True):
254 """Convert k-point density to Monkhorst-Pack grid size.
256 atoms: Atoms object
257 Contains unit cell and information about boundary conditions.
258 kptdensity: float
259 Required k-point density. Default value is 3.5 point per Ang^-1.
260 even: bool
261 Round up to even numbers.
262 """
264 recipcell = atoms.cell.reciprocal()
265 kpts = []
266 for i in range(3):
267 if atoms.pbc[i]:
268 k = 2 * pi * sqrt((recipcell[i] ** 2).sum()) * kptdensity
269 if even:
270 kpts.append(2 * int(np.ceil(k / 2)))
271 else:
272 kpts.append(int(np.ceil(k)))
273 else:
274 kpts.append(1)
275 return np.array(kpts)
278def kpts2mp(atoms, kpts, even=False):
279 if kpts is None:
280 return np.array([1, 1, 1])
281 if isinstance(kpts, (float, int)):
282 return kptdensity2monkhorstpack(atoms, kpts, even)
283 else:
284 return kpts
287def kpts2sizeandoffsets(
288 size=None, density=None, gamma=None, even=None, atoms=None
289):
290 """Helper function for selecting k-points.
292 Use either size or density.
294 size: 3 ints
295 Number of k-points.
296 density: float
297 K-point density in units of k-points per Ang^-1.
298 gamma: None or bool
299 Should the Gamma-point be included? Yes / no / don't care:
300 True / False / None.
301 even: None or bool
302 Should the number of k-points be even? Yes / no / don't care:
303 True / False / None.
304 atoms: Atoms object
305 Needed for calculating k-point density.
307 """
309 if size is not None and density is not None:
310 raise ValueError(
311 'Cannot specify k-point mesh size and ' 'density simultaneously'
312 )
313 elif density is not None and atoms is None:
314 raise ValueError(
315 'Cannot set k-points from "density" unless '
316 'Atoms are provided (need BZ dimensions).'
317 )
319 if size is None:
320 if density is None:
321 size = [1, 1, 1]
322 else:
323 size = kptdensity2monkhorstpack(atoms, density, None)
325 # Not using the rounding from kptdensity2monkhorstpack as it doesn't do
326 # rounding to odd numbers
327 if even is not None:
328 size = np.array(size)
329 remainder = size % 2
330 if even:
331 size += remainder
332 else: # Round up to odd numbers
333 size += 1 - remainder
335 offsets = [0, 0, 0]
336 if atoms is None:
337 pbc = [True, True, True]
338 else:
339 pbc = atoms.pbc
341 if gamma is not None:
342 for i, s in enumerate(size):
343 if pbc[i] and s % 2 != bool(gamma):
344 offsets[i] = 0.5 / s
346 return size, offsets
349@jsonable('kpoints')
350class KPoints:
351 def __init__(self, kpts=None):
352 if kpts is None:
353 kpts = np.zeros((1, 3))
354 self.kpts = kpts
356 def todict(self):
357 return vars(self)
360def kpts2kpts(kpts, atoms=None):
361 from ase.dft.kpoints import monkhorst_pack
363 if kpts is None:
364 return KPoints()
366 if hasattr(kpts, 'kpts'):
367 return kpts
369 if isinstance(kpts, dict):
370 if 'kpts' in kpts:
371 return KPoints(kpts['kpts'])
372 if 'path' in kpts:
373 cell = Cell.ascell(atoms.cell)
374 return cell.bandpath(pbc=atoms.pbc, **kpts)
375 size, offsets = kpts2sizeandoffsets(atoms=atoms, **kpts)
376 return KPoints(monkhorst_pack(size) + offsets)
378 if isinstance(kpts[0], int):
379 return KPoints(monkhorst_pack(kpts))
381 return KPoints(np.array(kpts))
384def kpts2ndarray(kpts, atoms=None):
385 """Convert kpts keyword to 2-d ndarray of scaled k-points."""
386 return kpts2kpts(kpts, atoms=atoms).kpts
389class EigenvalOccupationMixin:
390 """Define 'eigenvalues' and 'occupations' properties on class.
392 eigenvalues and occupations will be arrays of shape (spin, kpts, nbands).
394 Classes must implement the old-fashioned get_eigenvalues and
395 get_occupations methods."""
397 # We should maybe deprecate this and rely on the new
398 # Properties object for eigenvalues/occupations.
400 @property
401 def eigenvalues(self):
402 return self._propwrapper().eigenvalues
404 @property
405 def occupations(self):
406 return self._propwrapper().occupations
408 def _propwrapper(self):
409 from ase.calculator.singlepoint import OutputPropertyWrapper
411 return OutputPropertyWrapper(self)
414class Parameters(dict):
415 """Dictionary for parameters.
417 Special feature: If param is a Parameters instance, then param.xc
418 is a shorthand for param['xc'].
419 """
421 def __getattr__(self, key):
422 if key not in self:
423 return dict.__getattribute__(self, key)
424 return self[key]
426 def __setattr__(self, key, value):
427 self[key] = value
429 @classmethod
430 def read(cls, filename):
431 """Read parameters from file."""
432 # We use ast to evaluate literals, avoiding eval()
433 # for security reasons.
434 import ast
436 with open(filename) as fd:
437 txt = fd.read().strip()
438 assert txt.startswith('dict(')
439 assert txt.endswith(')')
440 txt = txt[5:-1]
442 # The tostring() representation "dict(...)" is not actually
443 # a literal, so we manually parse that along with the other
444 # formatting that we did manually:
445 dct = {}
446 for line in txt.splitlines():
447 key, val = line.split('=', 1)
448 key = key.strip()
449 val = val.strip()
450 if val[-1] == ',':
451 val = val[:-1]
452 dct[key] = ast.literal_eval(val)
454 parameters = cls(dct)
455 return parameters
457 def tostring(self):
458 keys = sorted(self)
459 return (
460 'dict('
461 + ',\n '.join(f'{key}={self[key]!r}' for key in keys)
462 + ')\n'
463 )
465 def write(self, filename):
466 Path(filename).write_text(self.tostring())
469class BaseCalculator(GetPropertiesMixin):
470 implemented_properties: List[str] = []
471 'Properties calculator can handle (energy, forces, ...)'
473 # Placeholder object for deprecated arguments. Let deprecated keywords
474 # default to _deprecated and then issue a warning if the user passed
475 # any other object (such as None).
476 _deprecated = object()
478 def __init__(self, parameters=None, use_cache=True):
479 if parameters is None:
480 parameters = {}
482 self.parameters = dict(parameters)
483 self.atoms = None
484 self.results = {}
485 self.use_cache = use_cache
487 def calculate_properties(self, atoms, properties):
488 """This method is experimental; currently for internal use."""
489 for name in properties:
490 if name not in all_outputs:
491 raise ValueError(f'No such property: {name}')
493 # We ignore system changes for now.
494 self.calculate(atoms, properties, system_changes=all_changes)
496 props = self.export_properties()
498 for name in properties:
499 if name not in props:
500 raise PropertyNotPresent(name)
501 return props
503 @abstractmethod
504 def calculate(self, atoms, properties, system_changes):
505 ...
507 def check_state(self, atoms, tol=1e-15):
508 """Check for any system changes since last calculation."""
509 if self.use_cache:
510 return compare_atoms(self.atoms, atoms, tol=tol)
511 else:
512 return all_changes
514 def get_property(self, name, atoms=None, allow_calculation=True):
515 if name not in self.implemented_properties:
516 raise PropertyNotImplementedError(
517 f'{name} property not implemented'
518 )
520 if atoms is None:
521 atoms = self.atoms
522 system_changes = []
523 else:
524 system_changes = self.check_state(atoms)
526 if system_changes:
527 self.atoms = None
528 self.results = {}
530 if name not in self.results:
531 if not allow_calculation:
532 return None
534 if self.use_cache:
535 self.atoms = atoms.copy()
537 self.calculate(atoms, [name], system_changes)
539 if name not in self.results:
540 # For some reason the calculator was not able to do what we want,
541 # and that is OK.
542 raise PropertyNotImplementedError(
543 '{} not present in this ' 'calculation'.format(name)
544 )
546 result = self.results[name]
547 if isinstance(result, np.ndarray):
548 result = result.copy()
549 return result
551 def calculation_required(self, atoms, properties):
552 assert not isinstance(properties, str)
553 system_changes = self.check_state(atoms)
554 if system_changes:
555 return True
556 for name in properties:
557 if name not in self.results:
558 return True
559 return False
561 def export_properties(self):
562 return Properties(self.results)
564 def _get_name(self) -> str: # child class can override this
565 return self.__class__.__name__.lower()
567 @property
568 def name(self) -> str:
569 return self._get_name()
572class Calculator(BaseCalculator):
573 """Base-class for all ASE calculators.
575 A calculator must raise PropertyNotImplementedError if asked for a
576 property that it can't calculate. So, if calculation of the
577 stress tensor has not been implemented, get_stress(atoms) should
578 raise PropertyNotImplementedError. This can be achieved simply by not
579 including the string 'stress' in the list implemented_properties
580 which is a class member. These are the names of the standard
581 properties: 'energy', 'forces', 'stress', 'dipole', 'charges',
582 'magmom' and 'magmoms'.
583 """
585 default_parameters: Dict[str, Any] = {}
586 'Default parameters'
588 ignored_changes: Set[str] = set()
589 'Properties of Atoms which we ignore for the purposes of cache '
590 'invalidation with check_state().'
592 discard_results_on_any_change = False
593 'Whether we purge the results following any change in the set() method. '
594 'Most (file I/O) calculators will probably want this.'
596 def __init__(
597 self,
598 restart=None,
599 ignore_bad_restart_file=BaseCalculator._deprecated,
600 label=None,
601 atoms=None,
602 directory='.',
603 **kwargs,
604 ):
605 """Basic calculator implementation.
607 restart: str
608 Prefix for restart file. May contain a directory. Default
609 is None: don't restart.
610 ignore_bad_restart_file: bool
611 Deprecated, please do not use.
612 Passing more than one positional argument to Calculator()
613 is deprecated and will stop working in the future.
614 Ignore broken or missing restart file. By default, it is an
615 error if the restart file is missing or broken.
616 directory: str or PurePath
617 Working directory in which to read and write files and
618 perform calculations.
619 label: str
620 Name used for all files. Not supported by all calculators.
621 May contain a directory, but please use the directory parameter
622 for that instead.
623 atoms: Atoms object
624 Optional Atoms object to which the calculator will be
625 attached. When restarting, atoms will get its positions and
626 unit-cell updated from file.
627 """
628 self.atoms = None # copy of atoms object from last calculation
629 self.results = {} # calculated properties (energy, forces, ...)
630 self.parameters = None # calculational parameters
631 self._directory = None # Initialize
633 if ignore_bad_restart_file is self._deprecated:
634 ignore_bad_restart_file = False
635 else:
636 warnings.warn(
637 FutureWarning(
638 'The keyword "ignore_bad_restart_file" is deprecated and '
639 'will be removed in a future version of ASE. Passing more '
640 'than one positional argument to Calculator is also '
641 'deprecated and will stop functioning in the future. '
642 'Please pass arguments by keyword (key=value) except '
643 'optionally the "restart" keyword.'
644 )
645 )
647 if restart is not None:
648 try:
649 self.read(restart) # read parameters, atoms and results
650 except ReadError:
651 if ignore_bad_restart_file:
652 self.reset()
653 else:
654 raise
656 self.directory = directory
657 self.prefix = None
658 if label is not None:
659 if self.directory == '.' and '/' in label:
660 # We specified directory in label, and nothing in the diretory
661 # key
662 self.label = label
663 elif '/' not in label:
664 # We specified our directory in the directory keyword
665 # or not at all
666 self.label = '/'.join((self.directory, label))
667 else:
668 raise ValueError(
669 'Directory redundantly specified though '
670 'directory="{}" and label="{}". '
671 'Please omit "/" in label.'.format(self.directory, label)
672 )
674 if self.parameters is None:
675 # Use default parameters if they were not read from file:
676 self.parameters = self.get_default_parameters()
678 if atoms is not None:
679 atoms.calc = self
680 if self.atoms is not None:
681 # Atoms were read from file. Update atoms:
682 if not (
683 equal(atoms.numbers, self.atoms.numbers)
684 and (atoms.pbc == self.atoms.pbc).all()
685 ):
686 raise CalculatorError('Atoms not compatible with file')
687 atoms.positions = self.atoms.positions
688 atoms.cell = self.atoms.cell
690 self.set(**kwargs)
692 if not hasattr(self, 'get_spin_polarized'):
693 self.get_spin_polarized = self._deprecated_get_spin_polarized
694 # XXX We are very naughty and do not call super constructor!
696 # For historical reasons we have a particular caching protocol.
697 # We disable the superclass' optional cache.
698 self.use_cache = False
700 @property
701 def directory(self) -> str:
702 return self._directory
704 @directory.setter
705 def directory(self, directory: Union[str, os.PathLike]):
706 self._directory = str(Path(directory)) # Normalize path.
708 @property
709 def label(self):
710 if self.directory == '.':
711 return self.prefix
713 # Generally, label ~ directory/prefix
714 #
715 # We use '/' rather than os.pathsep because
716 # 1) directory/prefix does not represent any actual path
717 # 2) We want the same string to work the same on all platforms
718 if self.prefix is None:
719 return self.directory + '/'
721 return f'{self.directory}/{self.prefix}'
723 @label.setter
724 def label(self, label):
725 if label is None:
726 self.directory = '.'
727 self.prefix = None
728 return
730 tokens = label.rsplit('/', 1)
731 if len(tokens) == 2:
732 directory, prefix = tokens
733 else:
734 assert len(tokens) == 1
735 directory = '.'
736 prefix = tokens[0]
737 if prefix == '':
738 prefix = None
739 self.directory = directory
740 self.prefix = prefix
742 def set_label(self, label):
743 """Set label and convert label to directory and prefix.
745 Examples:
747 * label='abc': (directory='.', prefix='abc')
748 * label='dir1/abc': (directory='dir1', prefix='abc')
749 * label=None: (directory='.', prefix=None)
750 """
751 self.label = label
753 def get_default_parameters(self):
754 return Parameters(copy.deepcopy(self.default_parameters))
756 def todict(self, skip_default=True):
757 defaults = self.get_default_parameters()
758 dct = {}
759 for key, value in self.parameters.items():
760 if hasattr(value, 'todict'):
761 value = value.todict()
762 if skip_default:
763 default = defaults.get(key, '_no_default_')
764 if default != '_no_default_' and equal(value, default):
765 continue
766 dct[key] = value
767 return dct
769 def reset(self):
770 """Clear all information from old calculation."""
772 self.atoms = None
773 self.results = {}
775 def read(self, label):
776 """Read atoms, parameters and calculated properties from output file.
778 Read result from self.label file. Raise ReadError if the file
779 is not there. If the file is corrupted or contains an error
780 message from the calculation, a ReadError should also be
781 raised. In case of succes, these attributes must set:
783 atoms: Atoms object
784 The state of the atoms from last calculation.
785 parameters: Parameters object
786 The parameter dictionary.
787 results: dict
788 Calculated properties like energy and forces.
790 The FileIOCalculator.read() method will typically read atoms
791 and parameters and get the results dict by calling the
792 read_results() method."""
794 self.set_label(label)
796 def get_atoms(self):
797 if self.atoms is None:
798 raise ValueError('Calculator has no atoms')
799 atoms = self.atoms.copy()
800 atoms.calc = self
801 return atoms
803 @classmethod
804 def read_atoms(cls, restart, **kwargs):
805 return cls(restart=restart, label=restart, **kwargs).get_atoms()
807 def set(self, **kwargs):
808 """Set parameters like set(key1=value1, key2=value2, ...).
810 A dictionary containing the parameters that have been changed
811 is returned.
813 Subclasses must implement a set() method that will look at the
814 chaneged parameters and decide if a call to reset() is needed.
815 If the changed parameters are harmless, like a change in
816 verbosity, then there is no need to call reset().
818 The special keyword 'parameters' can be used to read
819 parameters from a file."""
821 if 'parameters' in kwargs:
822 filename = kwargs.pop('parameters')
823 parameters = Parameters.read(filename)
824 parameters.update(kwargs)
825 kwargs = parameters
827 changed_parameters = {}
829 for key, value in kwargs.items():
830 oldvalue = self.parameters.get(key)
831 if key not in self.parameters or not equal(value, oldvalue):
832 changed_parameters[key] = value
833 self.parameters[key] = value
835 if self.discard_results_on_any_change and changed_parameters:
836 self.reset()
837 return changed_parameters
839 def check_state(self, atoms, tol=1e-15):
840 """Check for any system changes since last calculation."""
841 return compare_atoms(
842 self.atoms,
843 atoms,
844 tol=tol,
845 excluded_properties=set(self.ignored_changes),
846 )
848 def calculate(
849 self, atoms=None, properties=['energy'], system_changes=all_changes
850 ):
851 """Do the calculation.
853 properties: list of str
854 List of what needs to be calculated. Can be any combination
855 of 'energy', 'forces', 'stress', 'dipole', 'charges', 'magmom'
856 and 'magmoms'.
857 system_changes: list of str
858 List of what has changed since last calculation. Can be
859 any combination of these six: 'positions', 'numbers', 'cell',
860 'pbc', 'initial_charges' and 'initial_magmoms'.
862 Subclasses need to implement this, but can ignore properties
863 and system_changes if they want. Calculated properties should
864 be inserted into results dictionary like shown in this dummy
865 example::
867 self.results = {'energy': 0.0,
868 'forces': np.zeros((len(atoms), 3)),
869 'stress': np.zeros(6),
870 'dipole': np.zeros(3),
871 'charges': np.zeros(len(atoms)),
872 'magmom': 0.0,
873 'magmoms': np.zeros(len(atoms))}
875 The subclass implementation should first call this
876 implementation to set the atoms attribute and create any missing
877 directories.
878 """
879 if atoms is not None:
880 self.atoms = atoms.copy()
881 if not os.path.isdir(self._directory):
882 try:
883 os.makedirs(self._directory)
884 except FileExistsError as e:
885 # We can only end up here in case of a race condition if
886 # multiple Calculators are running concurrently *and* use the
887 # same _directory, which cannot be expected to work anyway.
888 msg = (
889 'Concurrent use of directory '
890 + self._directory
891 + 'by multiple Calculator instances detected. Please '
892 'use one directory per instance.'
893 )
894 raise RuntimeError(msg) from e
896 def calculate_numerical_forces(self, atoms, d=0.001):
897 """Calculate numerical forces using finite difference.
899 All atoms will be displaced by +d and -d in all directions."""
900 from ase.calculators.test import numeric_forces
902 return numeric_forces(atoms, d=d)
904 def calculate_numerical_stress(self, atoms, d=1e-6, voigt=True):
905 """Calculate numerical stress using finite difference."""
906 from ase.calculators.test import numeric_stress
908 return numeric_stress(atoms, d=d, voigt=voigt)
910 def _deprecated_get_spin_polarized(self):
911 msg = (
912 'This calculator does not implement get_spin_polarized(). '
913 'In the future, calc.get_spin_polarized() will work only on '
914 'calculator classes that explicitly implement this method or '
915 'inherit the method via specialized subclasses.'
916 )
917 warnings.warn(msg, FutureWarning)
918 return False
920 def band_structure(self):
921 """Create band-structure object for plotting."""
922 from ase.spectrum.band_structure import get_band_structure
924 # XXX This calculator is supposed to just have done a band structure
925 # calculation, but the calculator may not have the correct Fermi level
926 # if it updated the Fermi level after changing k-points.
927 # This will be a problem with some calculators (currently GPAW), and
928 # the user would have to override this by providing the Fermi level
929 # from the selfconsistent calculation.
930 return get_band_structure(calc=self)
933class OldShellProfile:
934 def __init__(self, name, command, prefix):
935 self.name = name
936 self.command = command
937 self.prefix = prefix
939 def execute(self, directory):
940 if self.command is None:
941 raise EnvironmentError(
942 'Please set ${} environment variable '.format(
943 'ASE_' + self.name.upper() + '_COMMAND'
944 )
945 + 'or supply the command keyword'
946 )
947 command = self.command
948 if 'PREFIX' in command:
949 command = command.replace('PREFIX', self.prefix)
951 try:
952 proc = subprocess.Popen(command, shell=True, cwd=directory)
953 except OSError as err:
954 # Actually this may never happen with shell=True, since
955 # probably the shell launches successfully. But we soon want
956 # to allow calling the subprocess directly, and then this
957 # distinction (failed to launch vs failed to run) is useful.
958 msg = f'Failed to execute "{command}"'
959 raise EnvironmentError(msg) from err
961 errorcode = proc.wait()
963 if errorcode:
964 path = os.path.abspath(directory)
965 msg = (
966 'Calculator "{}" failed with command "{}" failed in '
967 '{} with error code {}'.format(
968 self.name, command, path, errorcode
969 )
970 )
971 raise CalculationFailed(msg)
974class ArgvProfile:
975 def __init__(self, name, argv):
976 self.name = name
977 self.argv = argv
979 def execute(self, directory, stdout_name=None):
980 directory = Path(directory).resolve()
981 if stdout_name is None:
982 stdout_name = f'{self.name}.out'
983 stdout_path = directory / f'{stdout_name}.out'
984 try:
985 with open(stdout_path, 'w') as fd:
986 subprocess.run(self.argv, cwd=directory, check=True, stdout=fd)
987 except subprocess.CalledProcessError as err:
988 msg = (
989 f'Calculator {self.name} failed with args {self.argv} '
990 f'in directory {directory}'
991 )
992 raise CalculationFailed(msg) from err
993 return stdout_path
996class FileIOCalculator(Calculator):
997 """Base class for calculators that write/read input/output files."""
999 # command: Optional[str] = None
1000 # 'Command used to start calculation'
1002 # Fallback command when nothing else is specified.
1003 # There will be no fallback in the future; it must be explicitly
1004 # configured.
1005 _legacy_default_command: Optional[str] = None
1007 cfg = _cfg # Ensure easy access to config for subclasses
1009 def __init__(
1010 self,
1011 restart=None,
1012 ignore_bad_restart_file=Calculator._deprecated,
1013 label=None,
1014 atoms=None,
1015 command=None,
1016 profile=None,
1017 **kwargs,
1018 ):
1019 """File-IO calculator.
1021 command: str
1022 Command used to start calculation.
1023 """
1025 super().__init__(restart, ignore_bad_restart_file, label, atoms,
1026 **kwargs)
1028 if profile is None:
1029 profile = self._initialize_profile(command)
1030 self.profile = profile
1032 @property
1033 def command(self):
1034 # XXX deprecate me
1035 #
1036 # This is for calculators that invoke Popen directly on
1037 # self.command instead of lettung us (superclass) do it.
1038 return self.profile.command
1040 @command.setter
1041 def command(self, command):
1042 self.profile.command = command
1044 def _initialize_profile(self, command):
1045 if self.name in self.cfg.parser:
1046 section = self.cfg.parser[self.name]
1047 # XXX getargv() returns None if missing!
1048 return ArgvProfile(self.name, section.getargv('argv'))
1050 if command is None:
1051 name = 'ASE_' + self.name.upper() + '_COMMAND'
1052 command = self.cfg.get(name)
1054 if command is None:
1055 # XXX issue a FutureWarning if this causes the command
1056 # to no longer be None
1057 command = self._legacy_default_command
1059 if command is None:
1060 raise EnvironmentError(
1061 f'No configuration of {self.name}. '
1062 f'Missing section [{self.name}] in configuration')
1064 return OldShellProfile(self.name, command, self.prefix)
1066 def calculate(
1067 self, atoms=None, properties=['energy'], system_changes=all_changes
1068 ):
1069 Calculator.calculate(self, atoms, properties, system_changes)
1070 self.write_input(self.atoms, properties, system_changes)
1071 self.execute()
1072 self.read_results()
1074 def execute(self):
1075 self.profile.execute(self.directory)
1077 def write_input(self, atoms, properties=None, system_changes=None):
1078 """Write input file(s).
1080 Call this method first in subclasses so that directories are
1081 created automatically."""
1083 absdir = os.path.abspath(self.directory)
1084 if absdir != os.curdir and not os.path.isdir(self.directory):
1085 os.makedirs(self.directory)
1087 def read_results(self):
1088 """Read energy, forces, ... from output file(s)."""