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 'exciting': 'ExcitingGroundStateCalculator', 

139 'crystal': 'CRYSTAL', 

140 'ff': 'ForceField', 

141 'gamess_us': 'GAMESSUS', 

142 'gulp': 'GULP', 

143 'kim': 'KIM', 

144 'lammpsrun': 'LAMMPS', 

145 'lammpslib': 'LAMMPSlib', 

146 'lj': 'LennardJones', 

147 'mopac': 'MOPAC', 

148 'morse': 'MorsePotential', 

149 'nwchem': 'NWChem', 

150 'openmx': 'OpenMX', 

151 'orca': 'ORCA', 

152 'qchem': 'QChem', 

153 'tip3p': 'TIP3P', 

154 'tip4p': 'TIP4P'} 

155 

156 

157external_calculators = {} 

158 

159 

160def register_calculator_class(name, cls): 

161 """ Add the class into the database. """ 

162 assert name not in external_calculators 

163 external_calculators[name] = cls 

164 names.append(name) 

165 names.sort() 

166 

167 

168def get_calculator_class(name): 

169 """Return calculator class.""" 

170 if name == 'asap': 

171 from asap3 import EMT as Calculator 

172 elif name == 'gpaw': 

173 from gpaw import GPAW as Calculator 

174 elif name == 'hotbit': 

175 from hotbit import Calculator 

176 elif name == 'vasp2': 

177 from ase.calculators.vasp import Vasp2 as Calculator 

178 elif name == 'ace': 

179 from ase.calculators.acemolecule import ACE as Calculator 

180 elif name == 'Psi4': 

181 from ase.calculators.psi4 import Psi4 as Calculator 

182 elif name in external_calculators: 

183 Calculator = external_calculators[name] 

184 else: 

185 classname = special.get(name, name.title()) 

186 module = __import__('ase.calculators.' + name, {}, None, [classname]) 

187 Calculator = getattr(module, classname) 

188 return Calculator 

189 

190 

191def equal(a, b, tol=None, rtol=None, atol=None): 

192 """ndarray-enabled comparison function.""" 

193 # XXX Known bugs: 

194 # * Comparing cell objects (pbc not part of array representation) 

195 # * Infinite recursion for cyclic dicts 

196 # * Can of worms is open 

197 if tol is not None: 

198 msg = 'Use `equal(a, b, rtol=..., atol=...)` instead of `tol=...`' 

199 warnings.warn(msg, DeprecationWarning) 

200 assert rtol is None and atol is None, \ 

201 'Do not use deprecated `tol` with `atol` and/or `rtol`' 

202 rtol = tol 

203 atol = tol 

204 

205 a_is_dict = isinstance(a, dict) 

206 b_is_dict = isinstance(b, dict) 

207 if a_is_dict or b_is_dict: 

208 # Check that both a and b are dicts 

209 if not (a_is_dict and b_is_dict): 

210 return False 

211 if a.keys() != b.keys(): 

212 return False 

213 return all(equal(a[key], b[key], rtol=rtol, atol=atol) for key in a) 

214 

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

216 return False 

217 

218 if rtol is None and atol is None: 

219 return np.array_equal(a, b) 

220 

221 if rtol is None: 

222 rtol = 0 

223 if atol is None: 

224 atol = 0 

225 

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

227 

228 

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

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

231 

232 atoms: Atoms object 

233 Contains unit cell and information about boundary conditions. 

234 kptdensity: float 

235 Required k-point density. Default value is 3.5 point per Ang^-1. 

236 even: bool 

237 Round up to even numbers. 

238 """ 

239 

240 recipcell = atoms.cell.reciprocal() 

241 kpts = [] 

242 for i in range(3): 

243 if atoms.pbc[i]: 

244 k = 2 * pi * sqrt((recipcell[i]**2).sum()) * kptdensity 

245 if even: 

246 kpts.append(2 * int(np.ceil(k / 2))) 

247 else: 

248 kpts.append(int(np.ceil(k))) 

249 else: 

250 kpts.append(1) 

251 return np.array(kpts) 

252 

253 

254def kpts2mp(atoms, kpts, even=False): 

255 if kpts is None: 

256 return np.array([1, 1, 1]) 

257 if isinstance(kpts, (float, int)): 

258 return kptdensity2monkhorstpack(atoms, kpts, even) 

259 else: 

260 return kpts 

261 

262 

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

264 atoms=None): 

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

266 

267 Use either size or density. 

268 

269 size: 3 ints 

270 Number of k-points. 

271 density: float 

272 K-point density in units of k-points per Ang^-1. 

273 gamma: None or bool 

274 Should the Gamma-point be included? Yes / no / don't care: 

275 True / False / None. 

276 even: None or bool 

277 Should the number of k-points be even? Yes / no / don't care: 

278 True / False / None. 

279 atoms: Atoms object 

280 Needed for calculating k-point density. 

281 

282 """ 

283 

284 if size is not None and density is not None: 

285 raise ValueError('Cannot specify k-point mesh size and ' 

286 'density simultaneously') 

287 elif density is not None and atoms is None: 

288 raise ValueError('Cannot set k-points from "density" unless ' 

289 'Atoms are provided (need BZ dimensions).') 

290 

291 if size is None: 

292 if density is None: 

293 size = [1, 1, 1] 

294 else: 

295 size = kptdensity2monkhorstpack(atoms, density, None) 

296 

297 # Not using the rounding from kptdensity2monkhorstpack as it doesn't do 

298 # rounding to odd numbers 

299 if even is not None: 

300 size = np.array(size) 

301 remainder = size % 2 

302 if even: 

303 size += remainder 

304 else: # Round up to odd numbers 

305 size += (1 - remainder) 

306 

307 offsets = [0, 0, 0] 

308 if atoms is None: 

309 pbc = [True, True, True] 

310 else: 

311 pbc = atoms.pbc 

312 

313 if gamma is not None: 

314 for i, s in enumerate(size): 

315 if pbc[i] and s % 2 != bool(gamma): 

316 offsets[i] = 0.5 / s 

317 

318 return size, offsets 

319 

320 

321@jsonable('kpoints') 

322class KPoints: 

323 def __init__(self, kpts=None): 

324 if kpts is None: 

325 kpts = np.zeros((1, 3)) 

326 self.kpts = kpts 

327 

328 def todict(self): 

329 return vars(self) 

330 

331 

332def kpts2kpts(kpts, atoms=None): 

333 from ase.dft.kpoints import monkhorst_pack 

334 

335 if kpts is None: 

336 return KPoints() 

337 

338 if hasattr(kpts, 'kpts'): 

339 return kpts 

340 

341 if isinstance(kpts, dict): 

342 if 'kpts' in kpts: 

343 return KPoints(kpts['kpts']) 

344 if 'path' in kpts: 

345 cell = Cell.ascell(atoms.cell) 

346 return cell.bandpath(pbc=atoms.pbc, **kpts) 

347 size, offsets = kpts2sizeandoffsets(atoms=atoms, **kpts) 

348 return KPoints(monkhorst_pack(size) + offsets) 

349 

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

351 return KPoints(monkhorst_pack(kpts)) 

352 

353 return KPoints(np.array(kpts)) 

354 

355 

356def kpts2ndarray(kpts, atoms=None): 

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

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

359 

360 

361class EigenvalOccupationMixin: 

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

363 

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

365 

366 Classes must implement the old-fashioned get_eigenvalues and 

367 get_occupations methods.""" 

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

369 # Properties object for eigenvalues/occupations. 

370 

371 @property 

372 def eigenvalues(self): 

373 return self._propwrapper().eigenvalues 

374 

375 @property 

376 def occupations(self): 

377 return self._propwrapper().occupations 

378 

379 def _propwrapper(self): 

380 from ase.calculator.singlepoint import OutputPropertyWrapper 

381 return OutputPropertyWrapper(self) 

382 

383 

384class Parameters(dict): 

385 """Dictionary for parameters. 

386 

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

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

389 """ 

390 

391 def __getattr__(self, key): 

392 if key not in self: 

393 return dict.__getattribute__(self, key) 

394 return self[key] 

395 

396 def __setattr__(self, key, value): 

397 self[key] = value 

398 

399 @classmethod 

400 def read(cls, filename): 

401 """Read parameters from file.""" 

402 # We use ast to evaluate literals, avoiding eval() 

403 # for security reasons. 

404 import ast 

405 with open(filename) as fd: 

406 txt = fd.read().strip() 

407 assert txt.startswith('dict(') 

408 assert txt.endswith(')') 

409 txt = txt[5:-1] 

410 

411 # The tostring() representation "dict(...)" is not actually 

412 # a literal, so we manually parse that along with the other 

413 # formatting that we did manually: 

414 dct = {} 

415 for line in txt.splitlines(): 

416 key, val = line.split('=', 1) 

417 key = key.strip() 

418 val = val.strip() 

419 if val[-1] == ',': 

420 val = val[:-1] 

421 dct[key] = ast.literal_eval(val) 

422 

423 parameters = cls(dct) 

424 return parameters 

425 

426 def tostring(self): 

427 keys = sorted(self) 

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

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

430 

431 def write(self, filename): 

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

433 

434 

435class BaseCalculator(GetPropertiesMixin): 

436 implemented_properties: List[str] = [] 

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

438 

439 # Placeholder object for deprecated arguments. Let deprecated keywords 

440 # default to _deprecated and then issue a warning if the user passed 

441 # any other object (such as None). 

442 _deprecated = object() 

443 

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

445 if parameters is None: 

446 parameters = {} 

447 

448 self.parameters = dict(parameters) 

449 self.atoms = None 

450 self.results = {} 

451 self.use_cache = use_cache 

452 

453 def calculate_properties(self, atoms, properties): 

454 """This method is experimental; currently for internal use.""" 

455 for name in properties: 

456 if name not in all_outputs: 

457 raise ValueError(f'No such property: {name}') 

458 

459 # We ignore system changes for now. 

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

461 

462 props = self.export_properties() 

463 

464 for name in properties: 

465 if name not in props: 

466 raise PropertyNotPresent(name) 

467 return props 

468 

469 @abstractmethod 

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

471 ... 

472 

473 def check_state(self, atoms, tol=1e-15): 

474 """Check for any system changes since last calculation.""" 

475 if self.use_cache: 

476 return compare_atoms(self.atoms, atoms, tol=tol) 

477 else: 

478 return all_changes 

479 

480 def get_property(self, name, atoms=None, allow_calculation=True): 

481 if name not in self.implemented_properties: 

482 raise PropertyNotImplementedError('{} property not implemented' 

483 .format(name)) 

484 

485 if atoms is None: 

486 atoms = self.atoms 

487 system_changes = [] 

488 else: 

489 system_changes = self.check_state(atoms) 

490 

491 if system_changes: 

492 self.atoms = None 

493 self.results = {} 

494 

495 if name not in self.results: 

496 if not allow_calculation: 

497 return None 

498 

499 if self.use_cache: 

500 self.atoms = atoms.copy() 

501 

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

503 

504 if name not in self.results: 

505 # For some reason the calculator was not able to do what we want, 

506 # and that is OK. 

507 raise PropertyNotImplementedError('{} not present in this ' 

508 'calculation'.format(name)) 

509 

510 result = self.results[name] 

511 if isinstance(result, np.ndarray): 

512 result = result.copy() 

513 return result 

514 

515 def calculation_required(self, atoms, properties): 

516 assert not isinstance(properties, str) 

517 system_changes = self.check_state(atoms) 

518 if system_changes: 

519 return True 

520 for name in properties: 

521 if name not in self.results: 

522 return True 

523 return False 

524 

525 def export_properties(self): 

526 return Properties(self.results) 

527 

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

529 return self.__class__.__name__.lower() 

530 

531 @property 

532 def name(self) -> str: 

533 return self._get_name() 

534 

535 

536class Calculator(BaseCalculator): 

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

538 

539 A calculator must raise PropertyNotImplementedError if asked for a 

540 property that it can't calculate. So, if calculation of the 

541 stress tensor has not been implemented, get_stress(atoms) should 

542 raise PropertyNotImplementedError. This can be achieved simply by not 

543 including the string 'stress' in the list implemented_properties 

544 which is a class member. These are the names of the standard 

545 properties: 'energy', 'forces', 'stress', 'dipole', 'charges', 

546 'magmom' and 'magmoms'. 

547 """ 

548 

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

550 'Default parameters' 

551 

552 ignored_changes: Set[str] = set() 

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

554 'invalidation with check_state().' 

555 

556 discard_results_on_any_change = False 

557 'Whether we purge the results following any change in the set() method. ' 

558 'Most (file I/O) calculators will probably want this.' 

559 

560 def __init__(self, restart=None, 

561 ignore_bad_restart_file=BaseCalculator._deprecated, 

562 label=None, atoms=None, directory='.', 

563 **kwargs): 

564 """Basic calculator implementation. 

565 

566 restart: str 

567 Prefix for restart file. May contain a directory. Default 

568 is None: don't restart. 

569 ignore_bad_restart_file: bool 

570 Deprecated, please do not use. 

571 Passing more than one positional argument to Calculator() 

572 is deprecated and will stop working in the future. 

573 Ignore broken or missing restart file. By default, it is an 

574 error if the restart file is missing or broken. 

575 directory: str or PurePath 

576 Working directory in which to read and write files and 

577 perform calculations. 

578 label: str 

579 Name used for all files. Not supported by all calculators. 

580 May contain a directory, but please use the directory parameter 

581 for that instead. 

582 atoms: Atoms object 

583 Optional Atoms object to which the calculator will be 

584 attached. When restarting, atoms will get its positions and 

585 unit-cell updated from file. 

586 """ 

587 self.atoms = None # copy of atoms object from last calculation 

588 self.results = {} # calculated properties (energy, forces, ...) 

589 self.parameters = None # calculational parameters 

590 self._directory = None # Initialize 

591 

592 if ignore_bad_restart_file is self._deprecated: 

593 ignore_bad_restart_file = False 

594 else: 

595 warnings.warn(FutureWarning( 

596 'The keyword "ignore_bad_restart_file" is deprecated and ' 

597 'will be removed in a future version of ASE. Passing more ' 

598 'than one positional argument to Calculator is also ' 

599 'deprecated and will stop functioning in the future. ' 

600 'Please pass arguments by keyword (key=value) except ' 

601 'optionally the "restart" keyword.' 

602 )) 

603 

604 if restart is not None: 

605 try: 

606 self.read(restart) # read parameters, atoms and results 

607 except ReadError: 

608 if ignore_bad_restart_file: 

609 self.reset() 

610 else: 

611 raise 

612 

613 self.directory = directory 

614 self.prefix = None 

615 if label is not None: 

616 if self.directory == '.' and '/' in label: 

617 # We specified directory in label, and nothing in the diretory 

618 # key 

619 self.label = label 

620 elif '/' not in label: 

621 # We specified our directory in the directory keyword 

622 # or not at all 

623 self.label = '/'.join((self.directory, label)) 

624 else: 

625 raise ValueError('Directory redundantly specified though ' 

626 'directory="{}" and label="{}". ' 

627 'Please omit "/" in label.' 

628 .format(self.directory, label)) 

629 

630 if self.parameters is None: 

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

632 self.parameters = self.get_default_parameters() 

633 

634 if atoms is not None: 

635 atoms.calc = self 

636 if self.atoms is not None: 

637 # Atoms were read from file. Update atoms: 

638 if not (equal(atoms.numbers, self.atoms.numbers) and 

639 (atoms.pbc == self.atoms.pbc).all()): 

640 raise CalculatorError('Atoms not compatible with file') 

641 atoms.positions = self.atoms.positions 

642 atoms.cell = self.atoms.cell 

643 

644 self.set(**kwargs) 

645 

646 if not hasattr(self, 'get_spin_polarized'): 

647 self.get_spin_polarized = self._deprecated_get_spin_polarized 

648 # XXX We are very naughty and do not call super constructor! 

649 

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

651 # We disable the superclass' optional cache. 

652 self.use_cache = False 

653 

654 @property 

655 def directory(self) -> str: 

656 return self._directory 

657 

658 @directory.setter 

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

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

661 

662 @property 

663 def label(self): 

664 if self.directory == '.': 

665 return self.prefix 

666 

667 # Generally, label ~ directory/prefix 

668 # 

669 # We use '/' rather than os.pathsep because 

670 # 1) directory/prefix does not represent any actual path 

671 # 2) We want the same string to work the same on all platforms 

672 if self.prefix is None: 

673 return self.directory + '/' 

674 

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

676 

677 @label.setter 

678 def label(self, label): 

679 if label is None: 

680 self.directory = '.' 

681 self.prefix = None 

682 return 

683 

684 tokens = label.rsplit('/', 1) 

685 if len(tokens) == 2: 

686 directory, prefix = tokens 

687 else: 

688 assert len(tokens) == 1 

689 directory = '.' 

690 prefix = tokens[0] 

691 if prefix == '': 

692 prefix = None 

693 self.directory = directory 

694 self.prefix = prefix 

695 

696 def set_label(self, label): 

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

698 

699 Examples: 

700 

701 * label='abc': (directory='.', prefix='abc') 

702 * label='dir1/abc': (directory='dir1', prefix='abc') 

703 * label=None: (directory='.', prefix=None) 

704 """ 

705 self.label = label 

706 

707 def get_default_parameters(self): 

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

709 

710 def todict(self, skip_default=True): 

711 defaults = self.get_default_parameters() 

712 dct = {} 

713 for key, value in self.parameters.items(): 

714 if hasattr(value, 'todict'): 

715 value = value.todict() 

716 if skip_default: 

717 default = defaults.get(key, '_no_default_') 

718 if default != '_no_default_' and equal(value, default): 

719 continue 

720 dct[key] = value 

721 return dct 

722 

723 def reset(self): 

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

725 

726 self.atoms = None 

727 self.results = {} 

728 

729 def read(self, label): 

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

731 

732 Read result from self.label file. Raise ReadError if the file 

733 is not there. If the file is corrupted or contains an error 

734 message from the calculation, a ReadError should also be 

735 raised. In case of succes, these attributes must set: 

736 

737 atoms: Atoms object 

738 The state of the atoms from last calculation. 

739 parameters: Parameters object 

740 The parameter dictionary. 

741 results: dict 

742 Calculated properties like energy and forces. 

743 

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

745 and parameters and get the results dict by calling the 

746 read_results() method.""" 

747 

748 self.set_label(label) 

749 

750 def get_atoms(self): 

751 if self.atoms is None: 

752 raise ValueError('Calculator has no atoms') 

753 atoms = self.atoms.copy() 

754 atoms.calc = self 

755 return atoms 

756 

757 @classmethod 

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

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

760 

761 def set(self, **kwargs): 

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

763 

764 A dictionary containing the parameters that have been changed 

765 is returned. 

766 

767 Subclasses must implement a set() method that will look at the 

768 chaneged parameters and decide if a call to reset() is needed. 

769 If the changed parameters are harmless, like a change in 

770 verbosity, then there is no need to call reset(). 

771 

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

773 parameters from a file.""" 

774 

775 if 'parameters' in kwargs: 

776 filename = kwargs.pop('parameters') 

777 parameters = Parameters.read(filename) 

778 parameters.update(kwargs) 

779 kwargs = parameters 

780 

781 changed_parameters = {} 

782 

783 for key, value in kwargs.items(): 

784 oldvalue = self.parameters.get(key) 

785 if key not in self.parameters or not equal(value, oldvalue): 

786 changed_parameters[key] = value 

787 self.parameters[key] = value 

788 

789 if self.discard_results_on_any_change and changed_parameters: 

790 self.reset() 

791 return changed_parameters 

792 

793 def check_state(self, atoms, tol=1e-15): 

794 """Check for any system changes since last calculation.""" 

795 return compare_atoms(self.atoms, atoms, tol=tol, 

796 excluded_properties=set(self.ignored_changes)) 

797 

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

799 system_changes=all_changes): 

800 """Do the calculation. 

801 

802 properties: list of str 

803 List of what needs to be calculated. Can be any combination 

804 of 'energy', 'forces', 'stress', 'dipole', 'charges', 'magmom' 

805 and 'magmoms'. 

806 system_changes: list of str 

807 List of what has changed since last calculation. Can be 

808 any combination of these six: 'positions', 'numbers', 'cell', 

809 'pbc', 'initial_charges' and 'initial_magmoms'. 

810 

811 Subclasses need to implement this, but can ignore properties 

812 and system_changes if they want. Calculated properties should 

813 be inserted into results dictionary like shown in this dummy 

814 example:: 

815 

816 self.results = {'energy': 0.0, 

817 'forces': np.zeros((len(atoms), 3)), 

818 'stress': np.zeros(6), 

819 'dipole': np.zeros(3), 

820 'charges': np.zeros(len(atoms)), 

821 'magmom': 0.0, 

822 'magmoms': np.zeros(len(atoms))} 

823 

824 The subclass implementation should first call this 

825 implementation to set the atoms attribute and create any missing 

826 directories. 

827 """ 

828 if atoms is not None: 

829 self.atoms = atoms.copy() 

830 if not os.path.isdir(self._directory): 

831 try: 

832 os.makedirs(self._directory) 

833 except FileExistsError as e: 

834 # We can only end up here in case of a race condition if 

835 # multiple Calculators are running concurrently *and* use the 

836 # same _directory, which cannot be expected to work anyway. 

837 msg = ('Concurrent use of directory ' + self._directory + 

838 'by multiple Calculator instances detected. Please ' 

839 'use one directory per instance.') 

840 raise RuntimeError(msg) from e 

841 

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

843 """Calculate numerical forces using finite difference. 

844 

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

846 from ase.calculators.test import numeric_forces 

847 return numeric_forces(atoms, d=d) 

848 

849 def calculate_numerical_stress(self, atoms, d=1e-6, voigt=True): 

850 """Calculate numerical stress using finite difference.""" 

851 from ase.calculators.test import numeric_stress 

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

853 

854 def _deprecated_get_spin_polarized(self): 

855 msg = ('This calculator does not implement get_spin_polarized(). ' 

856 'In the future, calc.get_spin_polarized() will work only on ' 

857 'calculator classes that explicitly implement this method or ' 

858 'inherit the method via specialized subclasses.') 

859 warnings.warn(msg, FutureWarning) 

860 return False 

861 

862 def band_structure(self): 

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

864 from ase.spectrum.band_structure import get_band_structure 

865 # XXX This calculator is supposed to just have done a band structure 

866 # calculation, but the calculator may not have the correct Fermi level 

867 # if it updated the Fermi level after changing k-points. 

868 # This will be a problem with some calculators (currently GPAW), and 

869 # the user would have to override this by providing the Fermi level 

870 # from the selfconsistent calculation. 

871 return get_band_structure(calc=self) 

872 

873 

874class FileIOCalculator(Calculator): 

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

876 

877 command: Optional[str] = None 

878 'Command used to start calculation' 

879 

880 def __init__(self, restart=None, 

881 ignore_bad_restart_file=Calculator._deprecated, 

882 label=None, atoms=None, command=None, **kwargs): 

883 """File-IO calculator. 

884 

885 command: str 

886 Command used to start calculation. 

887 """ 

888 

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

890 atoms, **kwargs) 

891 

892 if command is not None: 

893 self.command = command 

894 else: 

895 name = 'ASE_' + self.name.upper() + '_COMMAND' 

896 self.command = os.environ.get(name, self.command) 

897 

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

899 system_changes=all_changes): 

900 Calculator.calculate(self, atoms, properties, system_changes) 

901 self.write_input(self.atoms, properties, system_changes) 

902 self.execute() 

903 self.read_results() 

904 

905 def execute(self): 

906 if self.command is None: 

907 raise CalculatorSetupError( 

908 'Please set ${} environment variable ' 

909 .format('ASE_' + self.name.upper() + '_COMMAND') + 

910 'or supply the command keyword') 

911 command = self.command 

912 if 'PREFIX' in command: 

913 command = command.replace('PREFIX', self.prefix) 

914 

915 try: 

916 proc = subprocess.Popen(command, shell=True, cwd=self.directory) 

917 except OSError as err: 

918 # Actually this may never happen with shell=True, since 

919 # probably the shell launches successfully. But we soon want 

920 # to allow calling the subprocess directly, and then this 

921 # distinction (failed to launch vs failed to run) is useful. 

922 msg = 'Failed to execute "{}"'.format(command) 

923 raise EnvironmentError(msg) from err 

924 

925 errorcode = proc.wait() 

926 

927 if errorcode: 

928 path = os.path.abspath(self.directory) 

929 msg = ('Calculator "{}" failed with command "{}" failed in ' 

930 '{} with error code {}'.format(self.name, command, 

931 path, errorcode)) 

932 raise CalculationFailed(msg) 

933 

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

935 """Write input file(s). 

936 

937 Call this method first in subclasses so that directories are 

938 created automatically.""" 

939 

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

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

942 os.makedirs(self.directory) 

943 

944 def read_results(self): 

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

946 pass