Coverage for /builds/debichem-team/python-ase/ase/calculators/calculator.py : 82.04%
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
5import pathlib
6from typing import Union, Optional, List, Set, Dict, Any
7import warnings
9import numpy as np
11from ase.cell import Cell
12from ase.outputs import Properties, all_outputs
13from ase.utils import jsonable
14from ase.calculators.abc import GetPropertiesMixin
17class CalculatorError(RuntimeError):
18 """Base class of error types related to ASE calculators."""
21class CalculatorSetupError(CalculatorError):
22 """Calculation cannot be performed with the given parameters.
24 Reasons to raise this errors are:
25 * The calculator is not properly configured
26 (missing executable, environment variables, ...)
27 * The given atoms object is not supported
28 * Calculator parameters are unsupported
30 Typically raised before a calculation."""
33class EnvironmentError(CalculatorSetupError):
34 """Raised if calculator is not properly set up with ASE.
36 May be missing an executable or environment variables."""
39class InputError(CalculatorSetupError):
40 """Raised if inputs given to the calculator were incorrect.
42 Bad input keywords or values, or missing pseudopotentials.
44 This may be raised before or during calculation, depending on
45 when the problem is detected."""
48class CalculationFailed(CalculatorError):
49 """Calculation failed unexpectedly.
51 Reasons to raise this error are:
52 * Calculation did not converge
53 * Calculation ran out of memory
54 * Segmentation fault or other abnormal termination
55 * Arithmetic trouble (singular matrices, NaN, ...)
57 Typically raised during calculation."""
60class SCFError(CalculationFailed):
61 """SCF loop did not converge."""
64class ReadError(CalculatorError):
65 """Unexpected irrecoverable error while reading calculation results."""
68class PropertyNotImplementedError(NotImplementedError):
69 """Raised if a calculator does not implement the requested property."""
72class PropertyNotPresent(CalculatorError):
73 """Requested property is missing.
75 Maybe it was never calculated, or for some reason was not extracted
76 with the rest of the results, without being a fatal ReadError."""
79def compare_atoms(atoms1, atoms2, tol=1e-15, excluded_properties=None):
80 """Check for system changes since last calculation. Properties in
81 ``excluded_properties`` are not checked."""
82 if atoms1 is None:
83 system_changes = all_changes[:]
84 else:
85 system_changes = []
87 properties_to_check = set(all_changes)
88 if excluded_properties:
89 properties_to_check -= set(excluded_properties)
91 # Check properties that aren't in Atoms.arrays but are attributes of
92 # Atoms objects
93 for prop in ['cell', 'pbc']:
94 if prop in properties_to_check:
95 properties_to_check.remove(prop)
96 if not equal(getattr(atoms1, prop), getattr(atoms2, prop),
97 atol=tol):
98 system_changes.append(prop)
100 arrays1 = set(atoms1.arrays)
101 arrays2 = set(atoms2.arrays)
103 # Add any properties that are only in atoms1.arrays or only in
104 # atoms2.arrays (and aren't excluded). Note that if, e.g. arrays1 has
105 # `initial_charges` which is merely zeros and arrays2 does not have
106 # this array, we'll still assume that the system has changed. However,
107 # this should only occur rarely.
108 system_changes += properties_to_check & (arrays1 ^ arrays2)
110 # Finally, check all of the non-excluded properties shared by the atoms
111 # arrays.
112 for prop in properties_to_check & arrays1 & arrays2:
113 if not equal(atoms1.arrays[prop], atoms2.arrays[prop], atol=tol):
114 system_changes.append(prop)
116 return system_changes
119all_properties = ['energy', 'forces', 'stress', 'stresses', 'dipole',
120 'charges', 'magmom', 'magmoms', 'free_energy', 'energies']
123all_changes = ['positions', 'numbers', 'cell', 'pbc',
124 'initial_charges', 'initial_magmoms']
127# Recognized names of calculators sorted alphabetically:
128names = ['abinit', 'ace', 'aims', 'amber', 'asap', 'castep', 'cp2k',
129 'crystal', 'demon', 'demonnano', 'dftb', 'dftd3', 'dmol', 'eam',
130 'elk', 'emt', 'espresso', 'exciting', 'ff', 'fleur', 'gamess_us',
131 'gaussian', 'gpaw', 'gromacs', 'gulp', 'hotbit', 'kim',
132 'lammpslib', 'lammpsrun', 'lj', 'mopac', 'morse', 'nwchem',
133 'octopus', 'onetep', 'openmx', 'orca', 'psi4', 'qchem', 'siesta',
134 'tip3p', 'tip4p', 'turbomole', 'vasp']
137special = {'cp2k': 'CP2K',
138 'demonnano': 'DemonNano',
139 'dftd3': 'DFTD3',
140 'dmol': 'DMol3',
141 'eam': 'EAM',
142 'elk': 'ELK',
143 'emt': 'EMT',
144 'crystal': 'CRYSTAL',
145 'ff': 'ForceField',
146 'fleur': 'FLEUR',
147 'gamess_us': 'GAMESSUS',
148 'gulp': 'GULP',
149 'kim': 'KIM',
150 'lammpsrun': 'LAMMPS',
151 'lammpslib': 'LAMMPSlib',
152 'lj': 'LennardJones',
153 'mopac': 'MOPAC',
154 'morse': 'MorsePotential',
155 'nwchem': 'NWChem',
156 'openmx': 'OpenMX',
157 'orca': 'ORCA',
158 'qchem': 'QChem',
159 'tip3p': 'TIP3P',
160 'tip4p': 'TIP4P'}
163external_calculators = {}
166def register_calculator_class(name, cls):
167 """ Add the class into the database. """
168 assert name not in external_calculators
169 external_calculators[name] = cls
170 names.append(name)
171 names.sort()
174def get_calculator_class(name):
175 """Return calculator class."""
176 if name == 'asap':
177 from asap3 import EMT as Calculator
178 elif name == 'gpaw':
179 from gpaw import GPAW as Calculator
180 elif name == 'hotbit':
181 from hotbit import Calculator
182 elif name == 'vasp2':
183 from ase.calculators.vasp import Vasp2 as Calculator
184 elif name == 'ace':
185 from ase.calculators.acemolecule import ACE as Calculator
186 elif name == 'Psi4':
187 from ase.calculators.psi4 import Psi4 as Calculator
188 elif name in external_calculators:
189 Calculator = external_calculators[name]
190 else:
191 classname = special.get(name, name.title())
192 module = __import__('ase.calculators.' + name, {}, None, [classname])
193 Calculator = getattr(module, classname)
194 return Calculator
197def equal(a, b, tol=None, rtol=None, atol=None):
198 """ndarray-enabled comparison function."""
199 # XXX Known bugs:
200 # * Comparing cell objects (pbc not part of array representation)
201 # * Infinite recursion for cyclic dicts
202 # * Can of worms is open
203 if tol is not None:
204 msg = 'Use `equal(a, b, rtol=..., atol=...)` instead of `tol=...`'
205 warnings.warn(msg, DeprecationWarning)
206 assert rtol is None and atol is None, \
207 'Do not use deprecated `tol` with `atol` and/or `rtol`'
208 rtol = tol
209 atol = tol
211 a_is_dict = isinstance(a, dict)
212 b_is_dict = isinstance(b, dict)
213 if a_is_dict or b_is_dict:
214 # Check that both a and b are dicts
215 if not (a_is_dict and b_is_dict):
216 return False
217 if a.keys() != b.keys():
218 return False
219 return all(equal(a[key], b[key], rtol=rtol, atol=atol) for key in a)
221 if np.shape(a) != np.shape(b):
222 return False
224 if rtol is None and atol is None:
225 return np.array_equal(a, b)
227 if rtol is None:
228 rtol = 0
229 if atol is None:
230 atol = 0
232 return np.allclose(a, b, rtol=rtol, atol=atol)
235def kptdensity2monkhorstpack(atoms, kptdensity=3.5, even=True):
236 """Convert k-point density to Monkhorst-Pack grid size.
238 atoms: Atoms object
239 Contains unit cell and information about boundary conditions.
240 kptdensity: float
241 Required k-point density. Default value is 3.5 point per Ang^-1.
242 even: bool
243 Round up to even numbers.
244 """
246 recipcell = atoms.cell.reciprocal()
247 kpts = []
248 for i in range(3):
249 if atoms.pbc[i]:
250 k = 2 * pi * sqrt((recipcell[i]**2).sum()) * kptdensity
251 if even:
252 kpts.append(2 * int(np.ceil(k / 2)))
253 else:
254 kpts.append(int(np.ceil(k)))
255 else:
256 kpts.append(1)
257 return np.array(kpts)
260def kpts2mp(atoms, kpts, even=False):
261 if kpts is None:
262 return np.array([1, 1, 1])
263 if isinstance(kpts, (float, int)):
264 return kptdensity2monkhorstpack(atoms, kpts, even)
265 else:
266 return kpts
269def kpts2sizeandoffsets(size=None, density=None, gamma=None, even=None,
270 atoms=None):
271 """Helper function for selecting k-points.
273 Use either size or density.
275 size: 3 ints
276 Number of k-points.
277 density: float
278 K-point density in units of k-points per Ang^-1.
279 gamma: None or bool
280 Should the Gamma-point be included? Yes / no / don't care:
281 True / False / None.
282 even: None or bool
283 Should the number of k-points be even? Yes / no / don't care:
284 True / False / None.
285 atoms: Atoms object
286 Needed for calculating k-point density.
288 """
290 if size is not None and density is not None:
291 raise ValueError('Cannot specify k-point mesh size and '
292 'density simultaneously')
293 elif density is not None and atoms is None:
294 raise ValueError('Cannot set k-points from "density" unless '
295 'Atoms are provided (need BZ dimensions).')
297 if size is None:
298 if density is None:
299 size = [1, 1, 1]
300 else:
301 size = kptdensity2monkhorstpack(atoms, density, None)
303 # Not using the rounding from kptdensity2monkhorstpack as it doesn't do
304 # rounding to odd numbers
305 if even is not None:
306 size = np.array(size)
307 remainder = size % 2
308 if even:
309 size += remainder
310 else: # Round up to odd numbers
311 size += (1 - remainder)
313 offsets = [0, 0, 0]
314 if atoms is None:
315 pbc = [True, True, True]
316 else:
317 pbc = atoms.pbc
319 if gamma is not None:
320 for i, s in enumerate(size):
321 if pbc[i] and s % 2 != bool(gamma):
322 offsets[i] = 0.5 / s
324 return size, offsets
327@jsonable('kpoints')
328class KPoints:
329 def __init__(self, kpts=None):
330 if kpts is None:
331 kpts = np.zeros((1, 3))
332 self.kpts = kpts
334 def todict(self):
335 return vars(self)
338def kpts2kpts(kpts, atoms=None):
339 from ase.dft.kpoints import monkhorst_pack
341 if kpts is None:
342 return KPoints()
344 if hasattr(kpts, 'kpts'):
345 return kpts
347 if isinstance(kpts, dict):
348 if 'kpts' in kpts:
349 return KPoints(kpts['kpts'])
350 if 'path' in kpts:
351 cell = Cell.ascell(atoms.cell)
352 return cell.bandpath(pbc=atoms.pbc, **kpts)
353 size, offsets = kpts2sizeandoffsets(atoms=atoms, **kpts)
354 return KPoints(monkhorst_pack(size) + offsets)
356 if isinstance(kpts[0], int):
357 return KPoints(monkhorst_pack(kpts))
359 return KPoints(np.array(kpts))
362def kpts2ndarray(kpts, atoms=None):
363 """Convert kpts keyword to 2-d ndarray of scaled k-points."""
364 return kpts2kpts(kpts, atoms=atoms).kpts
367class EigenvalOccupationMixin:
368 """Define 'eigenvalues' and 'occupations' properties on class.
370 eigenvalues and occupations will be arrays of shape (spin, kpts, nbands).
372 Classes must implement the old-fashioned get_eigenvalues and
373 get_occupations methods."""
375 @property
376 def eigenvalues(self):
377 return self.build_eig_occ_array(self.get_eigenvalues)
379 @property
380 def occupations(self):
381 return self.build_eig_occ_array(self.get_occupation_numbers)
383 def build_eig_occ_array(self, getter):
384 nspins = self.get_number_of_spins()
385 nkpts = len(self.get_ibz_k_points())
386 nbands = self.get_number_of_bands()
387 arr = np.zeros((nspins, nkpts, nbands))
388 for s in range(nspins):
389 for k in range(nkpts):
390 arr[s, k, :] = getter(spin=s, kpt=k)
391 return arr
394class Parameters(dict):
395 """Dictionary for parameters.
397 Special feature: If param is a Parameters instance, then param.xc
398 is a shorthand for param['xc'].
399 """
401 def __getattr__(self, key):
402 if key not in self:
403 return dict.__getattribute__(self, key)
404 return self[key]
406 def __setattr__(self, key, value):
407 self[key] = value
409 @classmethod
410 def read(cls, filename):
411 """Read parameters from file."""
412 # We use ast to evaluate literals, avoiding eval()
413 # for security reasons.
414 import ast
415 with open(filename) as fd:
416 txt = fd.read().strip()
417 assert txt.startswith('dict(')
418 assert txt.endswith(')')
419 txt = txt[5:-1]
421 # The tostring() representation "dict(...)" is not actually
422 # a literal, so we manually parse that along with the other
423 # formatting that we did manually:
424 dct = {}
425 for line in txt.splitlines():
426 key, val = line.split('=', 1)
427 key = key.strip()
428 val = val.strip()
429 if val[-1] == ',':
430 val = val[:-1]
431 dct[key] = ast.literal_eval(val)
433 parameters = cls(dct)
434 return parameters
436 def tostring(self):
437 keys = sorted(self)
438 return 'dict(' + ',\n '.join(
439 '{}={!r}'.format(key, self[key]) for key in keys) + ')\n'
441 def write(self, filename):
442 pathlib.Path(filename).write_text(self.tostring())
445class Calculator(GetPropertiesMixin):
446 """Base-class for all ASE calculators.
448 A calculator must raise PropertyNotImplementedError if asked for a
449 property that it can't calculate. So, if calculation of the
450 stress tensor has not been implemented, get_stress(atoms) should
451 raise PropertyNotImplementedError. This can be achieved simply by not
452 including the string 'stress' in the list implemented_properties
453 which is a class member. These are the names of the standard
454 properties: 'energy', 'forces', 'stress', 'dipole', 'charges',
455 'magmom' and 'magmoms'.
456 """
458 implemented_properties: List[str] = []
459 'Properties calculator can handle (energy, forces, ...)'
461 default_parameters: Dict[str, Any] = {}
462 'Default parameters'
464 ignored_changes: Set[str] = set()
465 'Properties of Atoms which we ignore for the purposes of cache '
466 'invalidation with check_state().'
468 discard_results_on_any_change = False
469 'Whether we purge the results following any change in the set() method. '
470 'Most (file I/O) calculators will probably want this.'
472 _deprecated = object()
474 def __init__(self, restart=None, ignore_bad_restart_file=_deprecated,
475 label=None, atoms=None, directory='.',
476 **kwargs):
477 """Basic calculator implementation.
479 restart: str
480 Prefix for restart file. May contain a directory. Default
481 is None: don't restart.
482 ignore_bad_restart_file: bool
483 Deprecated, please do not use.
484 Passing more than one positional argument to Calculator()
485 is deprecated and will stop working in the future.
486 Ignore broken or missing restart file. By default, it is an
487 error if the restart file is missing or broken.
488 directory: str or PurePath
489 Working directory in which to read and write files and
490 perform calculations.
491 label: str
492 Name used for all files. Not supported by all calculators.
493 May contain a directory, but please use the directory parameter
494 for that instead.
495 atoms: Atoms object
496 Optional Atoms object to which the calculator will be
497 attached. When restarting, atoms will get its positions and
498 unit-cell updated from file.
499 """
500 self.atoms = None # copy of atoms object from last calculation
501 self.results = {} # calculated properties (energy, forces, ...)
502 self.parameters = None # calculational parameters
503 self._directory = None # Initialize
505 if ignore_bad_restart_file is self._deprecated:
506 ignore_bad_restart_file = False
507 else:
508 warnings.warn(FutureWarning(
509 'The keyword "ignore_bad_restart_file" is deprecated and '
510 'will be removed in a future version of ASE. Passing more '
511 'than one positional argument to Calculator is also '
512 'deprecated and will stop functioning in the future. '
513 'Please pass arguments by keyword (key=value) except '
514 'optionally the "restart" keyword.'
515 ))
517 if restart is not None:
518 try:
519 self.read(restart) # read parameters, atoms and results
520 except ReadError:
521 if ignore_bad_restart_file:
522 self.reset()
523 else:
524 raise
526 self.directory = directory
527 self.prefix = None
528 if label is not None:
529 if self.directory == '.' and '/' in label:
530 # We specified directory in label, and nothing in the diretory key
531 self.label = label
532 elif '/' not in label:
533 # We specified our directory in the directory keyword
534 # or not at all
535 self.label = '/'.join((self.directory, label))
536 else:
537 raise ValueError('Directory redundantly specified though '
538 'directory="{}" and label="{}". '
539 'Please omit "/" in label.'
540 .format(self.directory, label))
542 if self.parameters is None:
543 # Use default parameters if they were not read from file:
544 self.parameters = self.get_default_parameters()
546 if atoms is not None:
547 atoms.calc = self
548 if self.atoms is not None:
549 # Atoms were read from file. Update atoms:
550 if not (equal(atoms.numbers, self.atoms.numbers) and
551 (atoms.pbc == self.atoms.pbc).all()):
552 raise CalculatorError('Atoms not compatible with file')
553 atoms.positions = self.atoms.positions
554 atoms.cell = self.atoms.cell
556 self.set(**kwargs)
558 if not hasattr(self, 'name'):
559 self.name = self.__class__.__name__.lower()
561 if not hasattr(self, 'get_spin_polarized'):
562 self.get_spin_polarized = self._deprecated_get_spin_polarized
564 @property
565 def directory(self) -> str:
566 return self._directory
568 @directory.setter
569 def directory(self, directory: Union[str, pathlib.PurePath]):
570 self._directory = str(pathlib.Path(directory)) # Normalize path.
572 @property
573 def label(self):
574 if self.directory == '.':
575 return self.prefix
577 # Generally, label ~ directory/prefix
578 #
579 # We use '/' rather than os.pathsep because
580 # 1) directory/prefix does not represent any actual path
581 # 2) We want the same string to work the same on all platforms
582 if self.prefix is None:
583 return self.directory + '/'
585 return '{}/{}'.format(self.directory, self.prefix)
587 @label.setter
588 def label(self, label):
589 if label is None:
590 self.directory = '.'
591 self.prefix = None
592 return
594 tokens = label.rsplit('/', 1)
595 if len(tokens) == 2:
596 directory, prefix = tokens
597 else:
598 assert len(tokens) == 1
599 directory = '.'
600 prefix = tokens[0]
601 if prefix == '':
602 prefix = None
603 self.directory = directory
604 self.prefix = prefix
606 def set_label(self, label):
607 """Set label and convert label to directory and prefix.
609 Examples:
611 * label='abc': (directory='.', prefix='abc')
612 * label='dir1/abc': (directory='dir1', prefix='abc')
613 * label=None: (directory='.', prefix=None)
614 """
615 self.label = label
617 def get_default_parameters(self):
618 return Parameters(copy.deepcopy(self.default_parameters))
620 def todict(self, skip_default=True):
621 defaults = self.get_default_parameters()
622 dct = {}
623 for key, value in self.parameters.items():
624 if hasattr(value, 'todict'):
625 value = value.todict()
626 if skip_default:
627 default = defaults.get(key, '_no_default_')
628 if default != '_no_default_' and equal(value, default):
629 continue
630 dct[key] = value
631 return dct
633 def reset(self):
634 """Clear all information from old calculation."""
636 self.atoms = None
637 self.results = {}
639 def read(self, label):
640 """Read atoms, parameters and calculated properties from output file.
642 Read result from self.label file. Raise ReadError if the file
643 is not there. If the file is corrupted or contains an error
644 message from the calculation, a ReadError should also be
645 raised. In case of succes, these attributes must set:
647 atoms: Atoms object
648 The state of the atoms from last calculation.
649 parameters: Parameters object
650 The parameter dictionary.
651 results: dict
652 Calculated properties like energy and forces.
654 The FileIOCalculator.read() method will typically read atoms
655 and parameters and get the results dict by calling the
656 read_results() method."""
658 self.set_label(label)
660 def get_atoms(self):
661 if self.atoms is None:
662 raise ValueError('Calculator has no atoms')
663 atoms = self.atoms.copy()
664 atoms.calc = self
665 return atoms
667 @classmethod
668 def read_atoms(cls, restart, **kwargs):
669 return cls(restart=restart, label=restart, **kwargs).get_atoms()
671 def set(self, **kwargs):
672 """Set parameters like set(key1=value1, key2=value2, ...).
674 A dictionary containing the parameters that have been changed
675 is returned.
677 Subclasses must implement a set() method that will look at the
678 chaneged parameters and decide if a call to reset() is needed.
679 If the changed parameters are harmless, like a change in
680 verbosity, then there is no need to call reset().
682 The special keyword 'parameters' can be used to read
683 parameters from a file."""
685 if 'parameters' in kwargs:
686 filename = kwargs.pop('parameters')
687 parameters = Parameters.read(filename)
688 parameters.update(kwargs)
689 kwargs = parameters
691 changed_parameters = {}
693 for key, value in kwargs.items():
694 oldvalue = self.parameters.get(key)
695 if key not in self.parameters or not equal(value, oldvalue):
696 changed_parameters[key] = value
697 self.parameters[key] = value
699 if self.discard_results_on_any_change and changed_parameters:
700 self.reset()
701 return changed_parameters
703 def check_state(self, atoms, tol=1e-15):
704 """Check for any system changes since last calculation."""
705 return compare_atoms(self.atoms, atoms, tol=tol,
706 excluded_properties=set(self.ignored_changes))
708 def get_potential_energy(self, atoms=None, force_consistent=False):
709 energy = self.get_property('energy', atoms)
710 if force_consistent:
711 if 'free_energy' not in self.results:
712 name = self.__class__.__name__
713 # XXX but we don't know why the energy is not there.
714 # We should raise PropertyNotPresent. Discuss
715 raise PropertyNotImplementedError(
716 'Force consistent/free energy ("free_energy") '
717 'not provided by {0} calculator'.format(name))
718 return self.results['free_energy']
719 else:
720 return energy
722 def get_property(self, name, atoms=None, allow_calculation=True):
723 if name not in self.implemented_properties:
724 raise PropertyNotImplementedError('{} property not implemented'
725 .format(name))
727 if atoms is None:
728 atoms = self.atoms
729 system_changes = []
730 else:
731 system_changes = self.check_state(atoms)
732 if system_changes:
733 self.reset()
734 if name not in self.results:
735 if not allow_calculation:
736 return None
737 self.calculate(atoms, [name], system_changes)
739 if name not in self.results:
740 # For some reason the calculator was not able to do what we want,
741 # and that is OK.
742 raise PropertyNotImplementedError('{} not present in this '
743 'calculation'.format(name))
745 result = self.results[name]
746 if isinstance(result, np.ndarray):
747 result = result.copy()
748 return result
750 def calculation_required(self, atoms, properties):
751 assert not isinstance(properties, str)
752 system_changes = self.check_state(atoms)
753 if system_changes:
754 return True
755 for name in properties:
756 if name not in self.results:
757 return True
758 return False
760 def calculate(self, atoms=None, properties=['energy'],
761 system_changes=all_changes):
762 """Do the calculation.
764 properties: list of str
765 List of what needs to be calculated. Can be any combination
766 of 'energy', 'forces', 'stress', 'dipole', 'charges', 'magmom'
767 and 'magmoms'.
768 system_changes: list of str
769 List of what has changed since last calculation. Can be
770 any combination of these six: 'positions', 'numbers', 'cell',
771 'pbc', 'initial_charges' and 'initial_magmoms'.
773 Subclasses need to implement this, but can ignore properties
774 and system_changes if they want. Calculated properties should
775 be inserted into results dictionary like shown in this dummy
776 example::
778 self.results = {'energy': 0.0,
779 'forces': np.zeros((len(atoms), 3)),
780 'stress': np.zeros(6),
781 'dipole': np.zeros(3),
782 'charges': np.zeros(len(atoms)),
783 'magmom': 0.0,
784 'magmoms': np.zeros(len(atoms))}
786 The subclass implementation should first call this
787 implementation to set the atoms attribute and create any missing
788 directories.
789 """
791 if atoms is not None:
792 self.atoms = atoms.copy()
793 if not os.path.isdir(self._directory):
794 os.makedirs(self._directory)
796 def calculate_numerical_forces(self, atoms, d=0.001):
797 """Calculate numerical forces using finite difference.
799 All atoms will be displaced by +d and -d in all directions."""
801 from ase.calculators.test import numeric_force
802 return np.array([[numeric_force(atoms, a, i, d)
803 for i in range(3)] for a in range(len(atoms))])
805 def calculate_numerical_stress(self, atoms, d=1e-6, voigt=True):
806 """Calculate numerical stress using finite difference."""
808 stress = np.zeros((3, 3), dtype=float)
810 cell = atoms.cell.copy()
811 V = atoms.get_volume()
812 for i in range(3):
813 x = np.eye(3)
814 x[i, i] += d
815 atoms.set_cell(np.dot(cell, x), scale_atoms=True)
816 eplus = atoms.get_potential_energy(force_consistent=True)
818 x[i, i] -= 2 * d
819 atoms.set_cell(np.dot(cell, x), scale_atoms=True)
820 eminus = atoms.get_potential_energy(force_consistent=True)
822 stress[i, i] = (eplus - eminus) / (2 * d * V)
823 x[i, i] += d
825 j = i - 2
826 x[i, j] = d
827 x[j, i] = d
828 atoms.set_cell(np.dot(cell, x), scale_atoms=True)
829 eplus = atoms.get_potential_energy(force_consistent=True)
831 x[i, j] = -d
832 x[j, i] = -d
833 atoms.set_cell(np.dot(cell, x), scale_atoms=True)
834 eminus = atoms.get_potential_energy(force_consistent=True)
836 stress[i, j] = (eplus - eminus) / (4 * d * V)
837 stress[j, i] = stress[i, j]
838 atoms.set_cell(cell, scale_atoms=True)
840 if voigt:
841 return stress.flat[[0, 4, 8, 5, 2, 1]]
842 else:
843 return stress
845 def _deprecated_get_spin_polarized(self):
846 msg = ('This calculator does not implement get_spin_polarized(). '
847 'In the future, calc.get_spin_polarized() will work only on '
848 'calculator classes that explicitly implement this method or '
849 'inherit the method via specialized subclasses.')
850 warnings.warn(msg, FutureWarning)
851 return False
853 def band_structure(self):
854 """Create band-structure object for plotting."""
855 from ase.spectrum.band_structure import get_band_structure
856 # XXX This calculator is supposed to just have done a band structure
857 # calculation, but the calculator may not have the correct Fermi level
858 # if it updated the Fermi level after changing k-points.
859 # This will be a problem with some calculators (currently GPAW), and
860 # the user would have to override this by providing the Fermi level
861 # from the selfconsistent calculation.
862 return get_band_structure(calc=self)
864 def calculate_properties(self, atoms, properties):
865 """This method is experimental; currently for internal use."""
866 for name in properties:
867 if name not in all_outputs:
868 raise ValueError(f'No such property: {name}')
870 # We ignore system changes for now.
871 self.calculate(atoms, properties, system_changes=all_changes)
873 props = self.export_properties()
875 for name in properties:
876 if name not in props:
877 raise PropertyNotPresent(name)
878 return props
880 def export_properties(self):
881 return Properties(self.results)
884class FileIOCalculator(Calculator):
885 """Base class for calculators that write/read input/output files."""
887 command: Optional[str] = None
888 'Command used to start calculation'
890 def __init__(self, restart=None,
891 ignore_bad_restart_file=Calculator._deprecated,
892 label=None, atoms=None, command=None, **kwargs):
893 """File-IO calculator.
895 command: str
896 Command used to start calculation.
897 """
899 Calculator.__init__(self, restart, ignore_bad_restart_file, label,
900 atoms, **kwargs)
902 if command is not None:
903 self.command = command
904 else:
905 name = 'ASE_' + self.name.upper() + '_COMMAND'
906 self.command = os.environ.get(name, self.command)
908 def calculate(self, atoms=None, properties=['energy'],
909 system_changes=all_changes):
910 Calculator.calculate(self, atoms, properties, system_changes)
911 self.write_input(self.atoms, properties, system_changes)
912 if self.command is None:
913 raise CalculatorSetupError(
914 'Please set ${} environment variable '
915 .format('ASE_' + self.name.upper() + '_COMMAND') +
916 'or supply the command keyword')
917 command = self.command
918 if 'PREFIX' in command:
919 command = command.replace('PREFIX', self.prefix)
921 try:
922 proc = subprocess.Popen(command, shell=True, cwd=self.directory)
923 except OSError as err:
924 # Actually this may never happen with shell=True, since
925 # probably the shell launches successfully. But we soon want
926 # to allow calling the subprocess directly, and then this
927 # distinction (failed to launch vs failed to run) is useful.
928 msg = 'Failed to execute "{}"'.format(command)
929 raise EnvironmentError(msg) from err
931 errorcode = proc.wait()
933 if errorcode:
934 path = os.path.abspath(self.directory)
935 msg = ('Calculator "{}" failed with command "{}" failed in '
936 '{} with error code {}'.format(self.name, command,
937 path, errorcode))
938 raise CalculationFailed(msg)
940 self.read_results()
942 def write_input(self, atoms, properties=None, system_changes=None):
943 """Write input file(s).
945 Call this method first in subclasses so that directories are
946 created automatically."""
948 absdir = os.path.abspath(self.directory)
949 if absdir != os.curdir and not os.path.isdir(self.directory):
950 os.makedirs(self.directory)
952 def read_results(self):
953 """Read energy, forces, ... from output file(s)."""
954 pass