Hide keyboard shortcuts

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 

9 

10import numpy as np 

11 

12from ase.cell import Cell 

13from ase.outputs import Properties, all_outputs 

14from ase.utils import jsonable 

15from ase.calculators.abc import GetPropertiesMixin 

16 

17from .names import names 

18 

19 

20class CalculatorError(RuntimeError): 

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

22 

23 

24class CalculatorSetupError(CalculatorError): 

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

26 

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 

32 

33 Typically raised before a calculation.""" 

34 

35 

36class EnvironmentError(CalculatorSetupError): 

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

38 

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

40 

41 

42class InputError(CalculatorSetupError): 

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

44 

45 Bad input keywords or values, or missing pseudopotentials. 

46 

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

48 when the problem is detected.""" 

49 

50 

51class CalculationFailed(CalculatorError): 

52 """Calculation failed unexpectedly. 

53 

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

59 

60 Typically raised during calculation.""" 

61 

62 

63class SCFError(CalculationFailed): 

64 """SCF loop did not converge.""" 

65 

66 

67class ReadError(CalculatorError): 

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

69 

70 

71class PropertyNotImplementedError(NotImplementedError): 

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

73 

74 

75class PropertyNotPresent(CalculatorError): 

76 """Requested property is missing. 

77 

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

80 

81 

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

89 

90 properties_to_check = set(all_changes) 

91 if excluded_properties: 

92 properties_to_check -= set(excluded_properties) 

93 

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) 

102 

103 arrays1 = set(atoms1.arrays) 

104 arrays2 = set(atoms2.arrays) 

105 

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) 

112 

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) 

118 

119 return system_changes 

120 

121 

122all_properties = ['energy', 'forces', 'stress', 'stresses', 'dipole', 

123 'charges', 'magmom', 'magmoms', 'free_energy', 'energies', 

124 'dielectric_tensor', 'born_effective_charges', 'polarization'] 

125 

126 

127all_changes = ['positions', 'numbers', 'cell', 'pbc', 

128 'initial_charges', 'initial_magmoms'] 

129 

130 

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

154 

155 

156external_calculators = {} 

157 

158 

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

165 

166 

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 

188 

189 

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 

203 

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) 

213 

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

215 return False 

216 

217 if rtol is None and atol is None: 

218 return np.array_equal(a, b) 

219 

220 if rtol is None: 

221 rtol = 0 

222 if atol is None: 

223 atol = 0 

224 

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

226 

227 

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

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

230 

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

238 

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) 

251 

252 

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 

260 

261 

262def kpts2sizeandoffsets(size=None, density=None, gamma=None, even=None, 

263 atoms=None): 

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

265 

266 Use either size or density. 

267 

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. 

280 

281 """ 

282 

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

289 

290 if size is None: 

291 if density is None: 

292 size = [1, 1, 1] 

293 else: 

294 size = kptdensity2monkhorstpack(atoms, density, None) 

295 

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) 

305 

306 offsets = [0, 0, 0] 

307 if atoms is None: 

308 pbc = [True, True, True] 

309 else: 

310 pbc = atoms.pbc 

311 

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 

316 

317 return size, offsets 

318 

319 

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 

326 

327 def todict(self): 

328 return vars(self) 

329 

330 

331def kpts2kpts(kpts, atoms=None): 

332 from ase.dft.kpoints import monkhorst_pack 

333 

334 if kpts is None: 

335 return KPoints() 

336 

337 if hasattr(kpts, 'kpts'): 

338 return kpts 

339 

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) 

348 

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

350 return KPoints(monkhorst_pack(kpts)) 

351 

352 return KPoints(np.array(kpts)) 

353 

354 

355def kpts2ndarray(kpts, atoms=None): 

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

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

358 

359 

360class EigenvalOccupationMixin: 

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

362 

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

364 

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. 

369 

370 @property 

371 def eigenvalues(self): 

372 return self._propwrapper().eigenvalues 

373 

374 @property 

375 def occupations(self): 

376 return self._propwrapper().occupations 

377 

378 def _propwrapper(self): 

379 from ase.calculator.singlepoint import OutputPropertyWrapper 

380 return OutputPropertyWrapper(self) 

381 

382 

383class Parameters(dict): 

384 """Dictionary for parameters. 

385 

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

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

388 """ 

389 

390 def __getattr__(self, key): 

391 if key not in self: 

392 return dict.__getattribute__(self, key) 

393 return self[key] 

394 

395 def __setattr__(self, key, value): 

396 self[key] = value 

397 

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] 

409 

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) 

421 

422 parameters = cls(dct) 

423 return parameters 

424 

425 def tostring(self): 

426 keys = sorted(self) 

427 return 'dict(' + ',\n '.join( 

428 '{}={!r}'.format(key, self[key]) for key in keys) + ')\n' 

429 

430 def write(self, filename): 

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

432 

433 

434class BaseCalculator(GetPropertiesMixin): 

435 implemented_properties: List[str] = [] 

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

437 

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

442 

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

444 if parameters is None: 

445 parameters = {} 

446 

447 self.parameters = dict(parameters) 

448 self.atoms = None 

449 self.results = {} 

450 self.use_cache = use_cache 

451 

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

457 

458 # We ignore system changes for now. 

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

460 

461 props = self.export_properties() 

462 

463 for name in properties: 

464 if name not in props: 

465 raise PropertyNotPresent(name) 

466 return props 

467 

468 @abstractmethod 

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

470 ... 

471 

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 

478 

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

483 

484 if atoms is None: 

485 atoms = self.atoms 

486 system_changes = [] 

487 else: 

488 system_changes = self.check_state(atoms) 

489 

490 if system_changes: 

491 self.atoms = None 

492 self.results = {} 

493 

494 if name not in self.results: 

495 if not allow_calculation: 

496 return None 

497 

498 if self.use_cache: 

499 self.atoms = atoms.copy() 

500 

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

502 

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

508 

509 result = self.results[name] 

510 if isinstance(result, np.ndarray): 

511 result = result.copy() 

512 return result 

513 

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 

523 

524 def export_properties(self): 

525 return Properties(self.results) 

526 

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

528 return self.__class__.__name__.lower() 

529 

530 @property 

531 def name(self) -> str: 

532 return self._get_name() 

533 

534 

535class Calculator(BaseCalculator): 

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

537 

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

547 

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

549 'Default parameters' 

550 

551 ignored_changes: Set[str] = set() 

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

553 'invalidation with check_state().' 

554 

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

558 

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. 

564 

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 

590 

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

602 

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 

611 

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

628 

629 if self.parameters is None: 

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

631 self.parameters = self.get_default_parameters() 

632 

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 

642 

643 self.set(**kwargs) 

644 

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! 

648 

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

650 # We disable the superclass' optional cache. 

651 self.use_cache = False 

652 

653 @property 

654 def directory(self) -> str: 

655 return self._directory 

656 

657 @directory.setter 

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

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

660 

661 @property 

662 def label(self): 

663 if self.directory == '.': 

664 return self.prefix 

665 

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

673 

674 return '{}/{}'.format(self.directory, self.prefix) 

675 

676 @label.setter 

677 def label(self, label): 

678 if label is None: 

679 self.directory = '.' 

680 self.prefix = None 

681 return 

682 

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 

694 

695 def set_label(self, label): 

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

697 

698 Examples: 

699 

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 

705 

706 def get_default_parameters(self): 

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

708 

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 

721 

722 def reset(self): 

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

724 

725 self.atoms = None 

726 self.results = {} 

727 

728 def read(self, label): 

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

730 

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: 

735 

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. 

742 

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

744 and parameters and get the results dict by calling the 

745 read_results() method.""" 

746 

747 self.set_label(label) 

748 

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 

755 

756 @classmethod 

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

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

759 

760 def set(self, **kwargs): 

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

762 

763 A dictionary containing the parameters that have been changed 

764 is returned. 

765 

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

770 

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

772 parameters from a file.""" 

773 

774 if 'parameters' in kwargs: 

775 filename = kwargs.pop('parameters') 

776 parameters = Parameters.read(filename) 

777 parameters.update(kwargs) 

778 kwargs = parameters 

779 

780 changed_parameters = {} 

781 

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 

787 

788 if self.discard_results_on_any_change and changed_parameters: 

789 self.reset() 

790 return changed_parameters 

791 

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

796 

797 def calculate(self, atoms=None, properties=['energy'], 

798 system_changes=all_changes): 

799 """Do the calculation. 

800 

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

809 

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

814 

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

822 

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 

840 

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

842 """Calculate numerical forces using finite difference. 

843 

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) 

847 

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) 

852 

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 

860 

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) 

871 

872 

873class FileIOCalculator(Calculator): 

874 """Base class for calculators that write/read input/output files.""" 

875 

876 command: Optional[str] = None 

877 'Command used to start calculation' 

878 

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. 

883 

884 command: str 

885 Command used to start calculation. 

886 """ 

887 

888 Calculator.__init__(self, restart, ignore_bad_restart_file, label, 

889 atoms, **kwargs) 

890 

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) 

896 

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

903 

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) 

913 

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 

923 

924 errorcode = proc.wait() 

925 

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) 

932 

933 def write_input(self, atoms, properties=None, system_changes=None): 

934 """Write input file(s). 

935 

936 Call this method first in subclasses so that directories are 

937 created automatically.""" 

938 

939 absdir = os.path.abspath(self.directory) 

940 if absdir != os.curdir and not os.path.isdir(self.directory): 

941 os.makedirs(self.directory) 

942 

943 def read_results(self): 

944 """Read energy, forces, ... from output file(s).""" 

945 pass