Coverage for /builds/ase/ase/ase/calculators/calculator.py : 81.82%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1import os
2import copy
3import subprocess
4from math import pi, sqrt
5from pathlib import Path
6from typing import Union, Optional, List, Set, Dict, Any
7import warnings
8from abc import abstractmethod
10import numpy as np
12from ase.cell import Cell
13from ase.outputs import Properties, all_outputs
14from ase.utils import jsonable
15from ase.calculators.abc import GetPropertiesMixin
17from .names import names
20class CalculatorError(RuntimeError):
21 """Base class of error types related to ASE calculators."""
24class CalculatorSetupError(CalculatorError):
25 """Calculation cannot be performed with the given parameters.
27 Reasons to raise this errors are:
28 * The calculator is not properly configured
29 (missing executable, environment variables, ...)
30 * The given atoms object is not supported
31 * Calculator parameters are unsupported
33 Typically raised before a calculation."""
36class EnvironmentError(CalculatorSetupError):
37 """Raised if calculator is not properly set up with ASE.
39 May be missing an executable or environment variables."""
42class InputError(CalculatorSetupError):
43 """Raised if inputs given to the calculator were incorrect.
45 Bad input keywords or values, or missing pseudopotentials.
47 This may be raised before or during calculation, depending on
48 when the problem is detected."""
51class CalculationFailed(CalculatorError):
52 """Calculation failed unexpectedly.
54 Reasons to raise this error are:
55 * Calculation did not converge
56 * Calculation ran out of memory
57 * Segmentation fault or other abnormal termination
58 * Arithmetic trouble (singular matrices, NaN, ...)
60 Typically raised during calculation."""
63class SCFError(CalculationFailed):
64 """SCF loop did not converge."""
67class ReadError(CalculatorError):
68 """Unexpected irrecoverable error while reading calculation results."""
71class PropertyNotImplementedError(NotImplementedError):
72 """Raised if a calculator does not implement the requested property."""
75class PropertyNotPresent(CalculatorError):
76 """Requested property is missing.
78 Maybe it was never calculated, or for some reason was not extracted
79 with the rest of the results, without being a fatal ReadError."""
82def compare_atoms(atoms1, atoms2, tol=1e-15, excluded_properties=None):
83 """Check for system changes since last calculation. Properties in
84 ``excluded_properties`` are not checked."""
85 if atoms1 is None:
86 system_changes = all_changes[:]
87 else:
88 system_changes = []
90 properties_to_check = set(all_changes)
91 if excluded_properties:
92 properties_to_check -= set(excluded_properties)
94 # Check properties that aren't in Atoms.arrays but are attributes of
95 # Atoms objects
96 for prop in ['cell', 'pbc']:
97 if prop in properties_to_check:
98 properties_to_check.remove(prop)
99 if not equal(getattr(atoms1, prop), getattr(atoms2, prop),
100 atol=tol):
101 system_changes.append(prop)
103 arrays1 = set(atoms1.arrays)
104 arrays2 = set(atoms2.arrays)
106 # Add any properties that are only in atoms1.arrays or only in
107 # atoms2.arrays (and aren't excluded). Note that if, e.g. arrays1 has
108 # `initial_charges` which is merely zeros and arrays2 does not have
109 # this array, we'll still assume that the system has changed. However,
110 # this should only occur rarely.
111 system_changes += properties_to_check & (arrays1 ^ arrays2)
113 # Finally, check all of the non-excluded properties shared by the atoms
114 # arrays.
115 for prop in properties_to_check & arrays1 & arrays2:
116 if not equal(atoms1.arrays[prop], atoms2.arrays[prop], atol=tol):
117 system_changes.append(prop)
119 return system_changes
122all_properties = ['energy', 'forces', 'stress', 'stresses', 'dipole',
123 'charges', 'magmom', 'magmoms', 'free_energy', 'energies',
124 'dielectric_tensor', 'born_effective_charges', 'polarization']
127all_changes = ['positions', 'numbers', 'cell', 'pbc',
128 'initial_charges', 'initial_magmoms']
131special = {'cp2k': 'CP2K',
132 'demonnano': 'DemonNano',
133 'dftd3': 'DFTD3',
134 'dmol': 'DMol3',
135 'eam': 'EAM',
136 'elk': 'ELK',
137 'emt': 'EMT',
138 'crystal': 'CRYSTAL',
139 'ff': 'ForceField',
140 'gamess_us': 'GAMESSUS',
141 'gulp': 'GULP',
142 'kim': 'KIM',
143 'lammpsrun': 'LAMMPS',
144 'lammpslib': 'LAMMPSlib',
145 'lj': 'LennardJones',
146 'mopac': 'MOPAC',
147 'morse': 'MorsePotential',
148 'nwchem': 'NWChem',
149 'openmx': 'OpenMX',
150 'orca': 'ORCA',
151 'qchem': 'QChem',
152 'tip3p': 'TIP3P',
153 'tip4p': 'TIP4P'}
156external_calculators = {}
159def register_calculator_class(name, cls):
160 """ Add the class into the database. """
161 assert name not in external_calculators
162 external_calculators[name] = cls
163 names.append(name)
164 names.sort()
167def get_calculator_class(name):
168 """Return calculator class."""
169 if name == 'asap':
170 from asap3 import EMT as Calculator
171 elif name == 'gpaw':
172 from gpaw import GPAW as Calculator
173 elif name == 'hotbit':
174 from hotbit import Calculator
175 elif name == 'vasp2':
176 from ase.calculators.vasp import Vasp2 as Calculator
177 elif name == 'ace':
178 from ase.calculators.acemolecule import ACE as Calculator
179 elif name == 'Psi4':
180 from ase.calculators.psi4 import Psi4 as Calculator
181 elif name in external_calculators:
182 Calculator = external_calculators[name]
183 else:
184 classname = special.get(name, name.title())
185 module = __import__('ase.calculators.' + name, {}, None, [classname])
186 Calculator = getattr(module, classname)
187 return Calculator
190def equal(a, b, tol=None, rtol=None, atol=None):
191 """ndarray-enabled comparison function."""
192 # XXX Known bugs:
193 # * Comparing cell objects (pbc not part of array representation)
194 # * Infinite recursion for cyclic dicts
195 # * Can of worms is open
196 if tol is not None:
197 msg = 'Use `equal(a, b, rtol=..., atol=...)` instead of `tol=...`'
198 warnings.warn(msg, DeprecationWarning)
199 assert rtol is None and atol is None, \
200 'Do not use deprecated `tol` with `atol` and/or `rtol`'
201 rtol = tol
202 atol = tol
204 a_is_dict = isinstance(a, dict)
205 b_is_dict = isinstance(b, dict)
206 if a_is_dict or b_is_dict:
207 # Check that both a and b are dicts
208 if not (a_is_dict and b_is_dict):
209 return False
210 if a.keys() != b.keys():
211 return False
212 return all(equal(a[key], b[key], rtol=rtol, atol=atol) for key in a)
214 if np.shape(a) != np.shape(b):
215 return False
217 if rtol is None and atol is None:
218 return np.array_equal(a, b)
220 if rtol is None:
221 rtol = 0
222 if atol is None:
223 atol = 0
225 return np.allclose(a, b, rtol=rtol, atol=atol)
228def kptdensity2monkhorstpack(atoms, kptdensity=3.5, even=True):
229 """Convert k-point density to Monkhorst-Pack grid size.
231 atoms: Atoms object
232 Contains unit cell and information about boundary conditions.
233 kptdensity: float
234 Required k-point density. Default value is 3.5 point per Ang^-1.
235 even: bool
236 Round up to even numbers.
237 """
239 recipcell = atoms.cell.reciprocal()
240 kpts = []
241 for i in range(3):
242 if atoms.pbc[i]:
243 k = 2 * pi * sqrt((recipcell[i]**2).sum()) * kptdensity
244 if even:
245 kpts.append(2 * int(np.ceil(k / 2)))
246 else:
247 kpts.append(int(np.ceil(k)))
248 else:
249 kpts.append(1)
250 return np.array(kpts)
253def kpts2mp(atoms, kpts, even=False):
254 if kpts is None:
255 return np.array([1, 1, 1])
256 if isinstance(kpts, (float, int)):
257 return kptdensity2monkhorstpack(atoms, kpts, even)
258 else:
259 return kpts
262def kpts2sizeandoffsets(size=None, density=None, gamma=None, even=None,
263 atoms=None):
264 """Helper function for selecting k-points.
266 Use either size or density.
268 size: 3 ints
269 Number of k-points.
270 density: float
271 K-point density in units of k-points per Ang^-1.
272 gamma: None or bool
273 Should the Gamma-point be included? Yes / no / don't care:
274 True / False / None.
275 even: None or bool
276 Should the number of k-points be even? Yes / no / don't care:
277 True / False / None.
278 atoms: Atoms object
279 Needed for calculating k-point density.
281 """
283 if size is not None and density is not None:
284 raise ValueError('Cannot specify k-point mesh size and '
285 'density simultaneously')
286 elif density is not None and atoms is None:
287 raise ValueError('Cannot set k-points from "density" unless '
288 'Atoms are provided (need BZ dimensions).')
290 if size is None:
291 if density is None:
292 size = [1, 1, 1]
293 else:
294 size = kptdensity2monkhorstpack(atoms, density, None)
296 # Not using the rounding from kptdensity2monkhorstpack as it doesn't do
297 # rounding to odd numbers
298 if even is not None:
299 size = np.array(size)
300 remainder = size % 2
301 if even:
302 size += remainder
303 else: # Round up to odd numbers
304 size += (1 - remainder)
306 offsets = [0, 0, 0]
307 if atoms is None:
308 pbc = [True, True, True]
309 else:
310 pbc = atoms.pbc
312 if gamma is not None:
313 for i, s in enumerate(size):
314 if pbc[i] and s % 2 != bool(gamma):
315 offsets[i] = 0.5 / s
317 return size, offsets
320@jsonable('kpoints')
321class KPoints:
322 def __init__(self, kpts=None):
323 if kpts is None:
324 kpts = np.zeros((1, 3))
325 self.kpts = kpts
327 def todict(self):
328 return vars(self)
331def kpts2kpts(kpts, atoms=None):
332 from ase.dft.kpoints import monkhorst_pack
334 if kpts is None:
335 return KPoints()
337 if hasattr(kpts, 'kpts'):
338 return kpts
340 if isinstance(kpts, dict):
341 if 'kpts' in kpts:
342 return KPoints(kpts['kpts'])
343 if 'path' in kpts:
344 cell = Cell.ascell(atoms.cell)
345 return cell.bandpath(pbc=atoms.pbc, **kpts)
346 size, offsets = kpts2sizeandoffsets(atoms=atoms, **kpts)
347 return KPoints(monkhorst_pack(size) + offsets)
349 if isinstance(kpts[0], int):
350 return KPoints(monkhorst_pack(kpts))
352 return KPoints(np.array(kpts))
355def kpts2ndarray(kpts, atoms=None):
356 """Convert kpts keyword to 2-d ndarray of scaled k-points."""
357 return kpts2kpts(kpts, atoms=atoms).kpts
360class EigenvalOccupationMixin:
361 """Define 'eigenvalues' and 'occupations' properties on class.
363 eigenvalues and occupations will be arrays of shape (spin, kpts, nbands).
365 Classes must implement the old-fashioned get_eigenvalues and
366 get_occupations methods."""
367 # We should maybe deprecate this and rely on the new
368 # Properties object for eigenvalues/occupations.
370 @property
371 def eigenvalues(self):
372 return self._propwrapper().eigenvalues
374 @property
375 def occupations(self):
376 return self._propwrapper().occupations
378 def _propwrapper(self):
379 from ase.calculator.singlepoint import OutputPropertyWrapper
380 return OutputPropertyWrapper(self)
383class Parameters(dict):
384 """Dictionary for parameters.
386 Special feature: If param is a Parameters instance, then param.xc
387 is a shorthand for param['xc'].
388 """
390 def __getattr__(self, key):
391 if key not in self:
392 return dict.__getattribute__(self, key)
393 return self[key]
395 def __setattr__(self, key, value):
396 self[key] = value
398 @classmethod
399 def read(cls, filename):
400 """Read parameters from file."""
401 # We use ast to evaluate literals, avoiding eval()
402 # for security reasons.
403 import ast
404 with open(filename) as fd:
405 txt = fd.read().strip()
406 assert txt.startswith('dict(')
407 assert txt.endswith(')')
408 txt = txt[5:-1]
410 # The tostring() representation "dict(...)" is not actually
411 # a literal, so we manually parse that along with the other
412 # formatting that we did manually:
413 dct = {}
414 for line in txt.splitlines():
415 key, val = line.split('=', 1)
416 key = key.strip()
417 val = val.strip()
418 if val[-1] == ',':
419 val = val[:-1]
420 dct[key] = ast.literal_eval(val)
422 parameters = cls(dct)
423 return parameters
425 def tostring(self):
426 keys = sorted(self)
427 return 'dict(' + ',\n '.join(
428 '{}={!r}'.format(key, self[key]) for key in keys) + ')\n'
430 def write(self, filename):
431 Path(filename).write_text(self.tostring())
434class BaseCalculator(GetPropertiesMixin):
435 implemented_properties: List[str] = []
436 'Properties calculator can handle (energy, forces, ...)'
438 # Placeholder object for deprecated arguments. Let deprecated keywords
439 # default to _deprecated and then issue a warning if the user passed
440 # any other object (such as None).
441 _deprecated = object()
443 def __init__(self, parameters=None, use_cache=True):
444 if parameters is None:
445 parameters = {}
447 self.parameters = dict(parameters)
448 self.atoms = None
449 self.results = {}
450 self.use_cache = use_cache
452 def calculate_properties(self, atoms, properties):
453 """This method is experimental; currently for internal use."""
454 for name in properties:
455 if name not in all_outputs:
456 raise ValueError(f'No such property: {name}')
458 # We ignore system changes for now.
459 self.calculate(atoms, properties, system_changes=all_changes)
461 props = self.export_properties()
463 for name in properties:
464 if name not in props:
465 raise PropertyNotPresent(name)
466 return props
468 @abstractmethod
469 def calculate(self, atoms, properties, system_changes):
470 ...
472 def check_state(self, atoms, tol=1e-15):
473 """Check for any system changes since last calculation."""
474 if self.use_cache:
475 return compare_atoms(self.atoms, atoms, tol=tol)
476 else:
477 return all_changes
479 def get_property(self, name, atoms=None, allow_calculation=True):
480 if name not in self.implemented_properties:
481 raise PropertyNotImplementedError('{} property not implemented'
482 .format(name))
484 if atoms is None:
485 atoms = self.atoms
486 system_changes = []
487 else:
488 system_changes = self.check_state(atoms)
490 if system_changes:
491 self.atoms = None
492 self.results = {}
494 if name not in self.results:
495 if not allow_calculation:
496 return None
498 if self.use_cache:
499 self.atoms = atoms.copy()
501 self.calculate(atoms, [name], system_changes)
503 if name not in self.results:
504 # For some reason the calculator was not able to do what we want,
505 # and that is OK.
506 raise PropertyNotImplementedError('{} not present in this '
507 'calculation'.format(name))
509 result = self.results[name]
510 if isinstance(result, np.ndarray):
511 result = result.copy()
512 return result
514 def calculation_required(self, atoms, properties):
515 assert not isinstance(properties, str)
516 system_changes = self.check_state(atoms)
517 if system_changes:
518 return True
519 for name in properties:
520 if name not in self.results:
521 return True
522 return False
524 def export_properties(self):
525 return Properties(self.results)
527 def _get_name(self) -> str: # child class can override this
528 return self.__class__.__name__.lower()
530 @property
531 def name(self) -> str:
532 return self._get_name()
535class Calculator(BaseCalculator):
536 """Base-class for all ASE calculators.
538 A calculator must raise PropertyNotImplementedError if asked for a
539 property that it can't calculate. So, if calculation of the
540 stress tensor has not been implemented, get_stress(atoms) should
541 raise PropertyNotImplementedError. This can be achieved simply by not
542 including the string 'stress' in the list implemented_properties
543 which is a class member. These are the names of the standard
544 properties: 'energy', 'forces', 'stress', 'dipole', 'charges',
545 'magmom' and 'magmoms'.
546 """
548 default_parameters: Dict[str, Any] = {}
549 'Default parameters'
551 ignored_changes: Set[str] = set()
552 'Properties of Atoms which we ignore for the purposes of cache '
553 'invalidation with check_state().'
555 discard_results_on_any_change = False
556 'Whether we purge the results following any change in the set() method. '
557 'Most (file I/O) calculators will probably want this.'
559 def __init__(self, restart=None,
560 ignore_bad_restart_file=BaseCalculator._deprecated,
561 label=None, atoms=None, directory='.',
562 **kwargs):
563 """Basic calculator implementation.
565 restart: str
566 Prefix for restart file. May contain a directory. Default
567 is None: don't restart.
568 ignore_bad_restart_file: bool
569 Deprecated, please do not use.
570 Passing more than one positional argument to Calculator()
571 is deprecated and will stop working in the future.
572 Ignore broken or missing restart file. By default, it is an
573 error if the restart file is missing or broken.
574 directory: str or PurePath
575 Working directory in which to read and write files and
576 perform calculations.
577 label: str
578 Name used for all files. Not supported by all calculators.
579 May contain a directory, but please use the directory parameter
580 for that instead.
581 atoms: Atoms object
582 Optional Atoms object to which the calculator will be
583 attached. When restarting, atoms will get its positions and
584 unit-cell updated from file.
585 """
586 self.atoms = None # copy of atoms object from last calculation
587 self.results = {} # calculated properties (energy, forces, ...)
588 self.parameters = None # calculational parameters
589 self._directory = None # Initialize
591 if ignore_bad_restart_file is self._deprecated:
592 ignore_bad_restart_file = False
593 else:
594 warnings.warn(FutureWarning(
595 'The keyword "ignore_bad_restart_file" is deprecated and '
596 'will be removed in a future version of ASE. Passing more '
597 'than one positional argument to Calculator is also '
598 'deprecated and will stop functioning in the future. '
599 'Please pass arguments by keyword (key=value) except '
600 'optionally the "restart" keyword.'
601 ))
603 if restart is not None:
604 try:
605 self.read(restart) # read parameters, atoms and results
606 except ReadError:
607 if ignore_bad_restart_file:
608 self.reset()
609 else:
610 raise
612 self.directory = directory
613 self.prefix = None
614 if label is not None:
615 if self.directory == '.' and '/' in label:
616 # We specified directory in label, and nothing in the diretory
617 # key
618 self.label = label
619 elif '/' not in label:
620 # We specified our directory in the directory keyword
621 # or not at all
622 self.label = '/'.join((self.directory, label))
623 else:
624 raise ValueError('Directory redundantly specified though '
625 'directory="{}" and label="{}". '
626 'Please omit "/" in label.'
627 .format(self.directory, label))
629 if self.parameters is None:
630 # Use default parameters if they were not read from file:
631 self.parameters = self.get_default_parameters()
633 if atoms is not None:
634 atoms.calc = self
635 if self.atoms is not None:
636 # Atoms were read from file. Update atoms:
637 if not (equal(atoms.numbers, self.atoms.numbers) and
638 (atoms.pbc == self.atoms.pbc).all()):
639 raise CalculatorError('Atoms not compatible with file')
640 atoms.positions = self.atoms.positions
641 atoms.cell = self.atoms.cell
643 self.set(**kwargs)
645 if not hasattr(self, 'get_spin_polarized'):
646 self.get_spin_polarized = self._deprecated_get_spin_polarized
647 # XXX We are very naughty and do not call super constructor!
649 # For historical reasons we have a particular caching protocol.
650 # We disable the superclass' optional cache.
651 self.use_cache = False
653 @property
654 def directory(self) -> str:
655 return self._directory
657 @directory.setter
658 def directory(self, directory: Union[str, os.PathLike]):
659 self._directory = str(Path(directory)) # Normalize path.
661 @property
662 def label(self):
663 if self.directory == '.':
664 return self.prefix
666 # Generally, label ~ directory/prefix
667 #
668 # We use '/' rather than os.pathsep because
669 # 1) directory/prefix does not represent any actual path
670 # 2) We want the same string to work the same on all platforms
671 if self.prefix is None:
672 return self.directory + '/'
674 return '{}/{}'.format(self.directory, self.prefix)
676 @label.setter
677 def label(self, label):
678 if label is None:
679 self.directory = '.'
680 self.prefix = None
681 return
683 tokens = label.rsplit('/', 1)
684 if len(tokens) == 2:
685 directory, prefix = tokens
686 else:
687 assert len(tokens) == 1
688 directory = '.'
689 prefix = tokens[0]
690 if prefix == '':
691 prefix = None
692 self.directory = directory
693 self.prefix = prefix
695 def set_label(self, label):
696 """Set label and convert label to directory and prefix.
698 Examples:
700 * label='abc': (directory='.', prefix='abc')
701 * label='dir1/abc': (directory='dir1', prefix='abc')
702 * label=None: (directory='.', prefix=None)
703 """
704 self.label = label
706 def get_default_parameters(self):
707 return Parameters(copy.deepcopy(self.default_parameters))
709 def todict(self, skip_default=True):
710 defaults = self.get_default_parameters()
711 dct = {}
712 for key, value in self.parameters.items():
713 if hasattr(value, 'todict'):
714 value = value.todict()
715 if skip_default:
716 default = defaults.get(key, '_no_default_')
717 if default != '_no_default_' and equal(value, default):
718 continue
719 dct[key] = value
720 return dct
722 def reset(self):
723 """Clear all information from old calculation."""
725 self.atoms = None
726 self.results = {}
728 def read(self, label):
729 """Read atoms, parameters and calculated properties from output file.
731 Read result from self.label file. Raise ReadError if the file
732 is not there. If the file is corrupted or contains an error
733 message from the calculation, a ReadError should also be
734 raised. In case of succes, these attributes must set:
736 atoms: Atoms object
737 The state of the atoms from last calculation.
738 parameters: Parameters object
739 The parameter dictionary.
740 results: dict
741 Calculated properties like energy and forces.
743 The FileIOCalculator.read() method will typically read atoms
744 and parameters and get the results dict by calling the
745 read_results() method."""
747 self.set_label(label)
749 def get_atoms(self):
750 if self.atoms is None:
751 raise ValueError('Calculator has no atoms')
752 atoms = self.atoms.copy()
753 atoms.calc = self
754 return atoms
756 @classmethod
757 def read_atoms(cls, restart, **kwargs):
758 return cls(restart=restart, label=restart, **kwargs).get_atoms()
760 def set(self, **kwargs):
761 """Set parameters like set(key1=value1, key2=value2, ...).
763 A dictionary containing the parameters that have been changed
764 is returned.
766 Subclasses must implement a set() method that will look at the
767 chaneged parameters and decide if a call to reset() is needed.
768 If the changed parameters are harmless, like a change in
769 verbosity, then there is no need to call reset().
771 The special keyword 'parameters' can be used to read
772 parameters from a file."""
774 if 'parameters' in kwargs:
775 filename = kwargs.pop('parameters')
776 parameters = Parameters.read(filename)
777 parameters.update(kwargs)
778 kwargs = parameters
780 changed_parameters = {}
782 for key, value in kwargs.items():
783 oldvalue = self.parameters.get(key)
784 if key not in self.parameters or not equal(value, oldvalue):
785 changed_parameters[key] = value
786 self.parameters[key] = value
788 if self.discard_results_on_any_change and changed_parameters:
789 self.reset()
790 return changed_parameters
792 def check_state(self, atoms, tol=1e-15):
793 """Check for any system changes since last calculation."""
794 return compare_atoms(self.atoms, atoms, tol=tol,
795 excluded_properties=set(self.ignored_changes))
797 def calculate(self, atoms=None, properties=['energy'],
798 system_changes=all_changes):
799 """Do the calculation.
801 properties: list of str
802 List of what needs to be calculated. Can be any combination
803 of 'energy', 'forces', 'stress', 'dipole', 'charges', 'magmom'
804 and 'magmoms'.
805 system_changes: list of str
806 List of what has changed since last calculation. Can be
807 any combination of these six: 'positions', 'numbers', 'cell',
808 'pbc', 'initial_charges' and 'initial_magmoms'.
810 Subclasses need to implement this, but can ignore properties
811 and system_changes if they want. Calculated properties should
812 be inserted into results dictionary like shown in this dummy
813 example::
815 self.results = {'energy': 0.0,
816 'forces': np.zeros((len(atoms), 3)),
817 'stress': np.zeros(6),
818 'dipole': np.zeros(3),
819 'charges': np.zeros(len(atoms)),
820 'magmom': 0.0,
821 'magmoms': np.zeros(len(atoms))}
823 The subclass implementation should first call this
824 implementation to set the atoms attribute and create any missing
825 directories.
826 """
827 if atoms is not None:
828 self.atoms = atoms.copy()
829 if not os.path.isdir(self._directory):
830 try:
831 os.makedirs(self._directory)
832 except FileExistsError as e:
833 # We can only end up here in case of a race condition if
834 # multiple Calculators are running concurrently *and* use the
835 # same _directory, which cannot be expected to work anyway.
836 msg = ('Concurrent use of directory ' + self._directory +
837 'by multiple Calculator instances detected. Please '
838 'use one directory per instance.')
839 raise RuntimeError(msg) from e
841 def calculate_numerical_forces(self, atoms, d=0.001):
842 """Calculate numerical forces using finite difference.
844 All atoms will be displaced by +d and -d in all directions."""
845 from ase.calculators.test import numeric_forces
846 return numeric_forces(atoms, d=d)
848 def calculate_numerical_stress(self, atoms, d=1e-6, voigt=True):
849 """Calculate numerical stress using finite difference."""
850 from ase.calculators.test import numeric_stress
851 return numeric_stress(atoms, d=d, voigt=voigt)
853 def _deprecated_get_spin_polarized(self):
854 msg = ('This calculator does not implement get_spin_polarized(). '
855 'In the future, calc.get_spin_polarized() will work only on '
856 'calculator classes that explicitly implement this method or '
857 'inherit the method via specialized subclasses.')
858 warnings.warn(msg, FutureWarning)
859 return False
861 def band_structure(self):
862 """Create band-structure object for plotting."""
863 from ase.spectrum.band_structure import get_band_structure
864 # XXX This calculator is supposed to just have done a band structure
865 # calculation, but the calculator may not have the correct Fermi level
866 # if it updated the Fermi level after changing k-points.
867 # This will be a problem with some calculators (currently GPAW), and
868 # the user would have to override this by providing the Fermi level
869 # from the selfconsistent calculation.
870 return get_band_structure(calc=self)
873class FileIOCalculator(Calculator):
874 """Base class for calculators that write/read input/output files."""
876 command: Optional[str] = None
877 'Command used to start calculation'
879 def __init__(self, restart=None,
880 ignore_bad_restart_file=Calculator._deprecated,
881 label=None, atoms=None, command=None, **kwargs):
882 """File-IO calculator.
884 command: str
885 Command used to start calculation.
886 """
888 Calculator.__init__(self, restart, ignore_bad_restart_file, label,
889 atoms, **kwargs)
891 if command is not None:
892 self.command = command
893 else:
894 name = 'ASE_' + self.name.upper() + '_COMMAND'
895 self.command = os.environ.get(name, self.command)
897 def calculate(self, atoms=None, properties=['energy'],
898 system_changes=all_changes):
899 Calculator.calculate(self, atoms, properties, system_changes)
900 self.write_input(self.atoms, properties, system_changes)
901 self.execute()
902 self.read_results()
904 def execute(self):
905 if self.command is None:
906 raise CalculatorSetupError(
907 'Please set ${} environment variable '
908 .format('ASE_' + self.name.upper() + '_COMMAND') +
909 'or supply the command keyword')
910 command = self.command
911 if 'PREFIX' in command:
912 command = command.replace('PREFIX', self.prefix)
914 try:
915 proc = subprocess.Popen(command, shell=True, cwd=self.directory)
916 except OSError as err:
917 # Actually this may never happen with shell=True, since
918 # probably the shell launches successfully. But we soon want
919 # to allow calling the subprocess directly, and then this
920 # distinction (failed to launch vs failed to run) is useful.
921 msg = 'Failed to execute "{}"'.format(command)
922 raise EnvironmentError(msg) from err
924 errorcode = proc.wait()
926 if errorcode:
927 path = os.path.abspath(self.directory)
928 msg = ('Calculator "{}" failed with command "{}" failed in '
929 '{} with error code {}'.format(self.name, command,
930 path, errorcode))
931 raise CalculationFailed(msg)
933 def write_input(self, atoms, properties=None, system_changes=None):
934 """Write input file(s).
936 Call this method first in subclasses so that directories are
937 created automatically."""
939 absdir = os.path.abspath(self.directory)
940 if absdir != os.curdir and not os.path.isdir(self.directory):
941 os.makedirs(self.directory)
943 def read_results(self):
944 """Read energy, forces, ... from output file(s)."""
945 pass