Coverage for /builds/hweiske/ase/ase/calculators/calculator.py: 84.60%

1import copy 

2import os 

3import subprocess 

4import warnings 

5from abc import abstractmethod 

6from dataclasses import dataclass 

7from math import pi, sqrt 

8from pathlib import Path 

9from typing import Any, Dict, List, Optional, Sequence, Set, Union 

10 

11import numpy as np 

12 

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 

18 

19from .names import names 

20 

21 

22class CalculatorError(RuntimeError): 

23 """Base class of error types related to ASE calculators.""" 

24 

25 

26class CalculatorSetupError(CalculatorError): 

27 """Calculation cannot be performed with the given parameters. 

28 

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 

34 

35 Typically raised before a calculation.""" 

36 

37 

38class EnvironmentError(CalculatorSetupError): 

39 """Raised if calculator is not properly set up with ASE. 

40 

41 May be missing an executable or environment variables.""" 

42 

43 

44class InputError(CalculatorSetupError): 

45 """Raised if inputs given to the calculator were incorrect. 

46 

47 Bad input keywords or values, or missing pseudopotentials. 

48 

49 This may be raised before or during calculation, depending on 

50 when the problem is detected.""" 

51 

52 

53class CalculationFailed(CalculatorError): 

54 """Calculation failed unexpectedly. 

55 

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, ...) 

61 

62 Typically raised during calculation.""" 

63 

64 

65class SCFError(CalculationFailed): 

66 """SCF loop did not converge.""" 

67 

68 

69class ReadError(CalculatorError): 

70 """Unexpected irrecoverable error while reading calculation results.""" 

71 

72 

73class PropertyNotImplementedError(NotImplementedError): 

74 """Raised if a calculator does not implement the requested property.""" 

75 

76 

77class PropertyNotPresent(CalculatorError): 

78 """Requested property is missing. 

79 

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.""" 

82 

83 

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 = [] 

91 

92 properties_to_check = set(all_changes) 

93 if excluded_properties: 

94 properties_to_check -= set(excluded_properties) 

95 

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) 

105 

106 arrays1 = set(atoms1.arrays) 

107 arrays2 = set(atoms2.arrays) 

108 

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) 

115 

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) 

121 

122 return system_changes 

123 

124 

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] 

140 

141 

142all_changes = [ 

143 'positions', 

144 'numbers', 

145 'cell', 

146 'pbc', 

147 'initial_charges', 

148 'initial_magmoms', 

149] 

150 

151 

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} 

178 

179 

180external_calculators = {} 

181 

182 

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() 

189 

190 

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 

212 

213 

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 

228 

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) 

238 

239 if np.shape(a) != np.shape(b): 

240 return False 

241 

242 if rtol is None and atol is None: 

243 return np.array_equal(a, b) 

244 

245 if rtol is None: 

246 rtol = 0 

247 if atol is None: 

248 atol = 0 

249 

250 return np.allclose(a, b, rtol=rtol, atol=atol) 

251 

252 

253def kptdensity2monkhorstpack(atoms, kptdensity=3.5, even=True): 

254 """Convert k-point density to Monkhorst-Pack grid size. 

255 

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 """ 

263 

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) 

276 

277 

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 

285 

286 

287def kpts2sizeandoffsets( 

288 size=None, density=None, gamma=None, even=None, atoms=None 

289): 

290 """Helper function for selecting k-points. 

291 

292 Use either size or density. 

293 

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. 

306 

307 """ 

308 

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 ) 

318 

319 if size is None: 

320 if density is None: 

321 size = [1, 1, 1] 

322 else: 

323 size = kptdensity2monkhorstpack(atoms, density, None) 

324 

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 

334 

335 offsets = [0, 0, 0] 

336 if atoms is None: 

337 pbc = [True, True, True] 

338 else: 

339 pbc = atoms.pbc 

340 

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 

345 

346 return size, offsets 

347 

348 

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 

355 

356 def todict(self): 

357 return vars(self) 

358 

359 

360def kpts2kpts(kpts, atoms=None): 

361 from ase.dft.kpoints import monkhorst_pack 

362 

363 if kpts is None: 

364 return KPoints() 

365 

366 if hasattr(kpts, 'kpts'): 

367 return kpts 

368 

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) 

377 

378 if isinstance(kpts[0], int): 

379 return KPoints(monkhorst_pack(kpts)) 

380 

381 return KPoints(np.array(kpts)) 

382 

383 

384def kpts2ndarray(kpts, atoms=None): 

385 """Convert kpts keyword to 2-d ndarray of scaled k-points.""" 

386 return kpts2kpts(kpts, atoms=atoms).kpts 

387 

388 

389class EigenvalOccupationMixin: 

390 """Define 'eigenvalues' and 'occupations' properties on class. 

391 

392 eigenvalues and occupations will be arrays of shape (spin, kpts, nbands). 

393 

394 Classes must implement the old-fashioned get_eigenvalues and 

395 get_occupations methods.""" 

396 

397 # We should maybe deprecate this and rely on the new 

398 # Properties object for eigenvalues/occupations. 

399 

400 @property 

401 def eigenvalues(self): 

402 return self._propwrapper().eigenvalues 

403 

404 @property 

405 def occupations(self): 

406 return self._propwrapper().occupations 

407 

408 def _propwrapper(self): 

409 from ase.calculator.singlepoint import OutputPropertyWrapper 

410 

411 return OutputPropertyWrapper(self) 

412 

413 

414class Parameters(dict): 

415 """Dictionary for parameters. 

416 

417 Special feature: If param is a Parameters instance, then param.xc 

418 is a shorthand for param['xc']. 

419 """ 

420 

421 def __getattr__(self, key): 

422 if key not in self: 

423 return dict.__getattribute__(self, key) 

424 return self[key] 

425 

426 def __setattr__(self, key, value): 

427 self[key] = value 

428 

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 

435 

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] 

441 

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) 

453 

454 parameters = cls(dct) 

455 return parameters 

456 

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 ) 

464 

465 def write(self, filename): 

466 Path(filename).write_text(self.tostring()) 

467 

468 

469class BaseCalculator(GetPropertiesMixin): 

470 implemented_properties: List[str] = [] 

471 'Properties calculator can handle (energy, forces, ...)' 

472 

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() 

477 

478 def __init__(self, parameters=None, use_cache=True): 

479 if parameters is None: 

480 parameters = {} 

481 

482 self.parameters = dict(parameters) 

483 self.atoms = None 

484 self.results = {} 

485 self.use_cache = use_cache 

486 

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}') 

492 

493 # We ignore system changes for now. 

494 self.calculate(atoms, properties, system_changes=all_changes) 

495 

496 props = self.export_properties() 

497 

498 for name in properties: 

499 if name not in props: 

500 raise PropertyNotPresent(name) 

501 return props 

502 

503 @abstractmethod 

504 def calculate(self, atoms, properties, system_changes): 

505 ... 

506 

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 

513 

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 ) 

519 

520 if atoms is None: 

521 atoms = self.atoms 

522 system_changes = [] 

523 else: 

524 system_changes = self.check_state(atoms) 

525 

526 if system_changes: 

527 self.atoms = None 

528 self.results = {} 

529 

530 if name not in self.results: 

531 if not allow_calculation: 

532 return None 

533 

534 if self.use_cache: 

535 self.atoms = atoms.copy() 

536 

537 self.calculate(atoms, [name], system_changes) 

538 

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 ) 

545 

546 result = self.results[name] 

547 if isinstance(result, np.ndarray): 

548 result = result.copy() 

549 return result 

550 

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 

560 

561 def export_properties(self): 

562 return Properties(self.results) 

563 

564 def _get_name(self) -> str: # child class can override this 

565 return self.__class__.__name__.lower() 

566 

567 @property 

568 def name(self) -> str: 

569 return self._get_name() 

570 

571 

572class Calculator(BaseCalculator): 

573 """Base-class for all ASE calculators. 

574 

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 """ 

584 

585 default_parameters: Dict[str, Any] = {} 

586 'Default parameters' 

587 

588 ignored_changes: Set[str] = set() 

589 'Properties of Atoms which we ignore for the purposes of cache ' 

590 'invalidation with check_state().' 

591 

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.' 

595 

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. 

606 

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 

632 

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 ) 

646 

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 

655 

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 ) 

673 

674 if self.parameters is None: 

675 # Use default parameters if they were not read from file: 

676 self.parameters = self.get_default_parameters() 

677 

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 

689 

690 self.set(**kwargs) 

691 

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! 

695 

696 # For historical reasons we have a particular caching protocol. 

697 # We disable the superclass' optional cache. 

698 self.use_cache = False 

699 

700 @property 

701 def directory(self) -> str: 

702 return self._directory 

703 

704 @directory.setter 

705 def directory(self, directory: Union[str, os.PathLike]): 

706 self._directory = str(Path(directory)) # Normalize path. 

707 

708 @property 

709 def label(self): 

710 if self.directory == '.': 

711 return self.prefix 

712 

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 + '/' 

720 

721 return f'{self.directory}/{self.prefix}' 

722 

723 @label.setter 

724 def label(self, label): 

725 if label is None: 

726 self.directory = '.' 

727 self.prefix = None 

728 return 

729 

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 

741 

742 def set_label(self, label): 

743 """Set label and convert label to directory and prefix. 

744 

745 Examples: 

746 

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 

752 

753 def get_default_parameters(self): 

754 return Parameters(copy.deepcopy(self.default_parameters)) 

755 

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 

768 

769 def reset(self): 

770 """Clear all information from old calculation.""" 

771 

772 self.atoms = None 

773 self.results = {} 

774 

775 def read(self, label): 

776 """Read atoms, parameters and calculated properties from output file. 

777 

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: 

782 

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. 

789 

790 The FileIOCalculator.read() method will typically read atoms 

791 and parameters and get the results dict by calling the 

792 read_results() method.""" 

793 

794 self.set_label(label) 

795 

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 

802 

803 @classmethod 

804 def read_atoms(cls, restart, **kwargs): 

805 return cls(restart=restart, label=restart, **kwargs).get_atoms() 

806 

807 def set(self, **kwargs): 

808 """Set parameters like set(key1=value1, key2=value2, ...). 

809 

810 A dictionary containing the parameters that have been changed 

811 is returned. 

812 

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(). 

817 

818 The special keyword 'parameters' can be used to read 

819 parameters from a file.""" 

820 

821 if 'parameters' in kwargs: 

822 filename = kwargs.pop('parameters') 

823 parameters = Parameters.read(filename) 

824 parameters.update(kwargs) 

825 kwargs = parameters 

826 

827 changed_parameters = {} 

828 

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 

834 

835 if self.discard_results_on_any_change and changed_parameters: 

836 self.reset() 

837 return changed_parameters 

838 

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 ) 

847 

848 def calculate( 

849 self, atoms=None, properties=['energy'], system_changes=all_changes 

850 ): 

851 """Do the calculation. 

852 

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'. 

861 

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:: 

866 

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))} 

874 

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 

895 

896 def calculate_numerical_forces(self, atoms, d=0.001): 

897 """Calculate numerical forces using finite difference. 

898 

899 All atoms will be displaced by +d and -d in all directions.""" 

900 from ase.calculators.test import numeric_forces 

901 

902 return numeric_forces(atoms, d=d) 

903 

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 

907 

908 return numeric_stress(atoms, d=d, voigt=voigt) 

909 

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 

919 

920 def band_structure(self): 

921 """Create band-structure object for plotting.""" 

922 from ase.spectrum.band_structure import get_band_structure 

923 

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