1import os 2import copy 3import subprocess 4from math import pi, sqrt 5import pathlib 6from typing import Union, Optional, List, Set, Dict, Any 7import warnings 8 9import numpy as np 10 11from ase.cell import Cell 12from ase.outputs import Properties, all_outputs 13from ase.utils import jsonable 14from ase.calculators.abc import GetPropertiesMixin 15 16 17class CalculatorError(RuntimeError): 18 """Base class of error types related to ASE calculators.""" 19 20 21class CalculatorSetupError(CalculatorError): 22 """Calculation cannot be performed with the given parameters. 23 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 29 30 Typically raised before a calculation.""" 31 32 33class EnvironmentError(CalculatorSetupError): 34 """Raised if calculator is not properly set up with ASE. 35 36 May be missing an executable or environment variables.""" 37 38 39class InputError(CalculatorSetupError): 40 """Raised if inputs given to the calculator were incorrect. 41 42 Bad input keywords or values, or missing pseudopotentials. 43 44 This may be raised before or during calculation, depending on 45 when the problem is detected.""" 46 47 48class CalculationFailed(CalculatorError): 49 """Calculation failed unexpectedly. 50 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, ...) 56 57 Typically raised during calculation.""" 58 59 60class SCFError(CalculationFailed): 61 """SCF loop did not converge.""" 62 63 64class ReadError(CalculatorError): 65 """Unexpected irrecoverable error while reading calculation results.""" 66 67 68class PropertyNotImplementedError(NotImplementedError): 69 """Raised if a calculator does not implement the requested property.""" 70 71 72class PropertyNotPresent(CalculatorError): 73 """Requested property is missing. 74 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.""" 77 78 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 = [] 86 87 properties_to_check = set(all_changes) 88 if excluded_properties: 89 properties_to_check -= set(excluded_properties) 90 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) 99 100 arrays1 = set(atoms1.arrays) 101 arrays2 = set(atoms2.arrays) 102 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) 109 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) 115 116 return system_changes 117 118 119all_properties = ['energy', 'forces', 'stress', 'stresses', 'dipole', 120 'charges', 'magmom', 'magmoms', 'free_energy', 'energies'] 121 122 123all_changes = ['positions', 'numbers', 'cell', 'pbc', 124 'initial_charges', 'initial_magmoms'] 125 126 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'] 135 136 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'} 161 162 163external_calculators = {} 164 165 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() 172 173 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 195 196 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 210 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) 220 221 if np.shape(a) != np.shape(b): 222 return False 223 224 if rtol is None and atol is None: 225 return np.array_equal(a, b) 226 227 if rtol is None: 228 rtol = 0 229 if atol is None: 230 atol = 0 231 232 return np.allclose(a, b, rtol=rtol, atol=atol) 233 234 235def kptdensity2monkhorstpack(atoms, kptdensity=3.5, even=True): 236 """Convert k-point density to Monkhorst-Pack grid size. 237 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 """ 245 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) 258 259 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 267 268 269def kpts2sizeandoffsets(size=None, density=None, gamma=None, even=None, 270 atoms=None): 271 """Helper function for selecting k-points. 272 273 Use either size or density. 274 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. 287 288 """ 289 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).') 296 297 if size is None: 298 if density is None: 299 size = [1, 1, 1] 300 else: 301 size = kptdensity2monkhorstpack(atoms, density, None) 302 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) 312 313 offsets = [0, 0, 0] 314 if atoms is None: 315 pbc = [True, True, True] 316 else: 317 pbc = atoms.pbc 318 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 323 324 return size, offsets 325 326 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 333 334 def todict(self): 335 return vars(self) 336 337 338def kpts2kpts(kpts, atoms=None): 339 from ase.dft.kpoints import monkhorst_pack 340 341 if kpts is None: 342 return KPoints() 343 344 if hasattr(kpts, 'kpts'): 345 return kpts 346 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) 355 356 if isinstance(kpts[0], int): 357 return KPoints(monkhorst_pack(kpts)) 358 359 return KPoints(np.array(kpts)) 360 361 362def kpts2ndarray(kpts, atoms=None): 363 """Convert kpts keyword to 2-d ndarray of scaled k-points.""" 364 return kpts2kpts(kpts, atoms=atoms).kpts 365 366 367class EigenvalOccupationMixin: 368 """Define 'eigenvalues' and 'occupations' properties on class. 369 370 eigenvalues and occupations will be arrays of shape (spin, kpts, nbands). 371 372 Classes must implement the old-fashioned get_eigenvalues and 373 get_occupations methods.""" 374 375 @property 376 def eigenvalues(self): 377 return self.build_eig_occ_array(self.get_eigenvalues) 378 379 @property 380 def occupations(self): 381 return self.build_eig_occ_array(self.get_occupation_numbers) 382 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 392 393 394class Parameters(dict): 395 """Dictionary for parameters. 396 397 Special feature: If param is a Parameters instance, then param.xc 398 is a shorthand for param['xc']. 399 """ 400 401 def __getattr__(self, key): 402 if key not in self: 403 return dict.__getattribute__(self, key) 404 return self[key] 405 406 def __setattr__(self, key, value): 407 self[key] = value 408 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] 420 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) 432 433 parameters = cls(dct) 434 return parameters 435 436 def tostring(self): 437 keys = sorted(self) 438 return 'dict(' + ',\n '.join( 439 '{}={!r}'.format(key, self[key]) for key in keys) + ')\n' 440 441 def write(self, filename): 442 pathlib.Path(filename).write_text(self.tostring()) 443 444 445class Calculator(GetPropertiesMixin): 446 """Base-class for all ASE calculators. 447 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 """ 457 458 implemented_properties: List[str] = [] 459 'Properties calculator can handle (energy, forces, ...)' 460 461 default_parameters: Dict[str, Any] = {} 462 'Default parameters' 463 464 ignored_changes: Set[str] = set() 465 'Properties of Atoms which we ignore for the purposes of cache ' 466 'invalidation with check_state().' 467 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.' 471 472 _deprecated = object() 473 474 def __init__(self, restart=None, ignore_bad_restart_file=_deprecated, 475 label=None, atoms=None, directory='.', 476 **kwargs): 477 """Basic calculator implementation. 478 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 504 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 )) 516 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 525 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)) 541 542 if self.parameters is None: 543 # Use default parameters if they were not read from file: 544 self.parameters = self.get_default_parameters() 545 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 555 556 self.set(**kwargs) 557 558 if not hasattr(self, 'name'): 559 self.name = self.__class__.__name__.lower() 560 561 if not hasattr(self, 'get_spin_polarized'): 562 self.get_spin_polarized = self._deprecated_get_spin_polarized 563 564 @property 565 def directory(self) -> str: 566 return self._directory 567 568 @directory.setter 569 def directory(self, directory: Union[str, pathlib.PurePath]): 570 self._directory = str(pathlib.Path(directory)) # Normalize path. 571 572 @property 573 def label(self): 574 if self.directory == '.': 575 return self.prefix 576 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 + '/' 584 585 return '{}/{}'.format(self.directory, self.prefix) 586 587 @label.setter 588 def label(self, label): 589 if label is None: 590 self.directory = '.' 591 self.prefix = None 592 return 593 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 605 606 def set_label(self, label): 607 """Set label and convert label to directory and prefix. 608 609 Examples: 610 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 616 617 def get_default_parameters(self): 618 return Parameters(copy.deepcopy(self.default_parameters)) 619 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 632 633 def reset(self): 634 """Clear all information from old calculation.""" 635 636 self.atoms = None 637 self.results = {} 638 639 def read(self, label): 640 """Read atoms, parameters and calculated properties from output file. 641 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: 646 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. 653 654 The FileIOCalculator.read() method will typically read atoms 655 and parameters and get the results dict by calling the 656 read_results() method.""" 657 658 self.set_label(label) 659 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 666 667 @classmethod 668 def read_atoms(cls, restart, **kwargs): 669 return cls(restart=restart, label=restart, **kwargs).get_atoms() 670 671 def set(self, **kwargs): 672 """Set parameters like set(key1=value1, key2=value2, ...). 673 674 A dictionary containing the parameters that have been changed 675 is returned. 676 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(). 681 682 The special keyword 'parameters' can be used to read 683 parameters from a file.""" 684 685 if 'parameters' in kwargs: 686 filename = kwargs.pop('parameters') 687 parameters = Parameters.read(filename) 688 parameters.update(kwargs) 689 kwargs = parameters 690 691 changed_parameters = {} 692 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 698 699 if self.discard_results_on_any_change and changed_parameters: 700 self.reset() 701 return changed_parameters 702 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)) 707 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 721 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)) 726 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) 738 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)) 744 745 result = self.results[name] 746 if isinstance(result, np.ndarray): 747 result = result.copy() 748 return result 749 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 759 760 def calculate(self, atoms=None, properties=['energy'], 761 system_changes=all_changes): 762 """Do the calculation. 763 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'. 772 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:: 777 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))} 785 786 The subclass implementation should first call this 787 implementation to set the atoms attribute and create any missing 788 directories. 789 """ 790 791 if atoms is not None: 792 self.atoms = atoms.copy() 793 if not os.path.isdir(self._directory): 794 os.makedirs(self._directory) 795 796 def calculate_numerical_forces(self, atoms, d=0.001): 797 """Calculate numerical forces using finite difference. 798 799 All atoms will be displaced by +d and -d in all directions.""" 800 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))]) 804 805 def calculate_numerical_stress(self, atoms, d=1e-6, voigt=True): 806 """Calculate numerical stress using finite difference.""" 807 808 stress = np.zeros((3, 3), dtype=float) 809 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) 817 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) 821 822 stress[i, i] = (eplus - eminus) / (2 * d * V) 823 x[i, i] += d 824 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) 830 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) 835 836 stress[i, j] = (eplus - eminus) / (4 * d * V) 837 stress[j, i] = stress[i, j] 838 atoms.set_cell(cell, scale_atoms=True) 839 840 if voigt: 841 return stress.flat[[0, 4, 8, 5, 2, 1]] 842 else: 843 return stress 844 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 852 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) 863 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}') 869 870 # We ignore system changes for now. 871 self.calculate(atoms, properties, system_changes=all_changes) 872 873 props = self.export_properties() 874 875 for name in properties: 876 if name not in props: 877 raise PropertyNotPresent(name) 878 return props 879 880 def export_properties(self): 881 return Properties(self.results) 882 883 884class FileIOCalculator(Calculator): 885 """Base class for calculators that write/read input/output files.""" 886 887 command: Optional[str] = None 888 'Command used to start calculation' 889 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. 894 895 command: str 896 Command used to start calculation. 897 """ 898 899 Calculator.__init__(self, restart, ignore_bad_restart_file, label, 900 atoms, **kwargs) 901 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) 907 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) 920 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 930 931 errorcode = proc.wait() 932 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) 939 940 self.read_results() 941 942 def write_input(self, atoms, properties=None, system_changes=None): 943 """Write input file(s). 944 945 Call this method first in subclasses so that directories are 946 created automatically.""" 947 948 absdir = os.path.abspath(self.directory) 949 if absdir != os.curdir and not os.path.isdir(self.directory): 950 os.makedirs(self.directory) 951 952 def read_results(self): 953 """Read energy, forces, ... from output file(s).""" 954 pass 955