1# Copyright (C) 2008 CSC - Scientific Computing Ltd. 2"""This module defines an ASE interface to VASP. 3 4Developed on the basis of modules by Jussi Enkovaara and John 5Kitchin. The path of the directory containing the pseudopotential 6directories (potpaw,potpaw_GGA, potpaw_PBE, ...) should be set 7by the environmental flag $VASP_PP_PATH. 8 9The user should also set the environmental flag $VASP_SCRIPT pointing 10to a python script looking something like:: 11 12 import os 13 exitcode = os.system('vasp') 14 15Alternatively, user can set the environmental flag $VASP_COMMAND pointing 16to the command use the launch vasp e.g. 'vasp' or 'mpirun -n 16 vasp' 17 18http://cms.mpi.univie.ac.at/vasp/ 19""" 20 21import os 22import sys 23import re 24import numpy as np 25import subprocess 26from contextlib import contextmanager 27from pathlib import Path 28from warnings import warn 29from typing import Dict, Any 30from xml.etree import ElementTree 31 32import ase 33from ase.io import read, jsonio 34from ase.utils import PurePath 35from ase.calculators import calculator 36from ase.calculators.calculator import Calculator 37from ase.calculators.singlepoint import SinglePointDFTCalculator 38from ase.calculators.vasp.create_input import GenerateVaspInput 39 40 41class Vasp(GenerateVaspInput, Calculator): # type: ignore 42 """ASE interface for the Vienna Ab initio Simulation Package (VASP), 43 with the Calculator interface. 44 45 Parameters: 46 47 atoms: object 48 Attach an atoms object to the calculator. 49 50 label: str 51 Prefix for the output file, and sets the working directory. 52 Default is 'vasp'. 53 54 directory: str 55 Set the working directory. Is prepended to ``label``. 56 57 restart: str or bool 58 Sets a label for the directory to load files from. 59 if :code:`restart=True`, the working directory from 60 ``directory`` is used. 61 62 txt: bool, None, str or writable object 63 - If txt is None, output stream will be supressed 64 65 - If txt is '-' the output will be sent through stdout 66 67 - If txt is a string a file will be opened,\ 68 and the output will be sent to that file. 69 70 - Finally, txt can also be a an output stream,\ 71 which has a 'write' attribute. 72 73 Default is 'vasp.out' 74 75 - Examples: 76 77 >>> Vasp(label='mylabel', txt='vasp.out') # Redirect stdout 78 >>> Vasp(txt='myfile.txt') # Redirect stdout 79 >>> Vasp(txt='-') # Print vasp output to stdout 80 >>> Vasp(txt=None) # Suppress txt output 81 82 command: str 83 Custom instructions on how to execute VASP. Has priority over 84 environment variables. 85 """ 86 name = 'vasp' 87 ase_objtype = 'vasp_calculator' # For JSON storage 88 89 # Environment commands 90 env_commands = ('ASE_VASP_COMMAND', 'VASP_COMMAND', 'VASP_SCRIPT') 91 92 implemented_properties = [ 93 'energy', 'free_energy', 'forces', 'dipole', 'fermi', 'stress', 94 'magmom', 'magmoms' 95 ] 96 97 # Can be used later to set some ASE defaults 98 default_parameters: Dict[str, Any] = {} 99 100 def __init__(self, 101 atoms=None, 102 restart=None, 103 directory='.', 104 label='vasp', 105 ignore_bad_restart_file=Calculator._deprecated, 106 command=None, 107 txt='vasp.out', 108 **kwargs): 109 110 self._atoms = None 111 self.results = {} 112 113 # Initialize parameter dictionaries 114 GenerateVaspInput.__init__(self) 115 self._store_param_state() # Initialize an empty parameter state 116 117 # Store calculator from vasprun.xml here - None => uninitialized 118 self._xml_calc = None 119 120 # Set directory and label 121 self.directory = directory 122 if '/' in label: 123 warn(('Specifying directory in "label" is deprecated, ' 124 'use "directory" instead.'), np.VisibleDeprecationWarning) 125 if self.directory != '.': 126 raise ValueError('Directory redundantly specified though ' 127 'directory="{}" and label="{}". ' 128 'Please omit "/" in label.'.format( 129 self.directory, label)) 130 self.label = label 131 else: 132 self.prefix = label # The label should only contain the prefix 133 134 if isinstance(restart, bool): 135 if restart is True: 136 restart = self.label 137 else: 138 restart = None 139 140 Calculator.__init__( 141 self, 142 restart=restart, 143 ignore_bad_restart_file=ignore_bad_restart_file, 144 # We already, manually, created the label 145 label=self.label, 146 atoms=atoms, 147 **kwargs) 148 149 self.command = command 150 151 self._txt = None 152 self.txt = txt # Set the output txt stream 153 self.version = None 154 155 # XXX: This seems to break restarting, unless we return first. 156 # Do we really still need to enfore this? 157 158 # # If no XC combination, GGA functional or POTCAR type is specified, 159 # # default to PW91. This is mostly chosen for backwards compatibility. 160 # if kwargs.get('xc', None): 161 # pass 162 # elif not (kwargs.get('gga', None) or kwargs.get('pp', None)): 163 # self.input_params.update({'xc': 'PW91'}) 164 # # A null value of xc is permitted; custom recipes can be 165 # # used by explicitly setting the pseudopotential set and 166 # # INCAR keys 167 # else: 168 # self.input_params.update({'xc': None}) 169 170 def make_command(self, command=None): 171 """Return command if one is passed, otherwise try to find 172 ASE_VASP_COMMAND, VASP_COMMAND or VASP_SCRIPT. 173 If none are set, a CalculatorSetupError is raised""" 174 if command: 175 cmd = command 176 else: 177 # Search for the environment commands 178 for env in self.env_commands: 179 if env in os.environ: 180 cmd = os.environ[env].replace('PREFIX', self.prefix) 181 if env == 'VASP_SCRIPT': 182 # Make the system python exe run $VASP_SCRIPT 183 exe = sys.executable 184 cmd = ' '.join([exe, cmd]) 185 break 186 else: 187 msg = ('Please set either command in calculator' 188 ' or one of the following environment ' 189 'variables (prioritized as follows): {}').format( 190 ', '.join(self.env_commands)) 191 raise calculator.CalculatorSetupError(msg) 192 return cmd 193 194 def set(self, **kwargs): 195 """Override the set function, to test for changes in the 196 Vasp Calculator, then call the create_input.set() 197 on remaining inputs for VASP specific keys. 198 199 Allows for setting ``label``, ``directory`` and ``txt`` 200 without resetting the results in the calculator. 201 """ 202 changed_parameters = {} 203 204 if 'label' in kwargs: 205 self.label = kwargs.pop('label') 206 207 if 'directory' in kwargs: 208 # str() call to deal with pathlib objects 209 self.directory = str(kwargs.pop('directory')) 210 211 if 'txt' in kwargs: 212 self.txt = kwargs.pop('txt') 213 214 if 'atoms' in kwargs: 215 atoms = kwargs.pop('atoms') 216 self.atoms = atoms # Resets results 217 218 if 'command' in kwargs: 219 self.command = kwargs.pop('command') 220 221 changed_parameters.update(Calculator.set(self, **kwargs)) 222 223 # We might at some point add more to changed parameters, or use it 224 if changed_parameters: 225 self.clear_results() # We don't want to clear atoms 226 if kwargs: 227 # If we make any changes to Vasp input, we always reset 228 GenerateVaspInput.set(self, **kwargs) 229 self.results.clear() 230 231 def reset(self): 232 self.atoms = None 233 self.clear_results() 234 235 def clear_results(self): 236 self.results.clear() 237 self._xml_calc = None 238 239 @contextmanager 240 def _txt_outstream(self): 241 """Custom function for opening a text output stream. Uses self.txt 242 to determine the output stream, and accepts a string or an open 243 writable object. 244 If a string is used, a new stream is opened, and automatically closes 245 the new stream again when exiting. 246 247 Examples: 248 # Pass a string 249 calc.txt = 'vasp.out' 250 with calc.txt_outstream() as out: 251 calc.run(out=out) # Redirects the stdout to 'vasp.out' 252 253 # Use an existing stream 254 mystream = open('vasp.out', 'w') 255 calc.txt = mystream 256 with calc.txt_outstream() as out: 257 calc.run(out=out) 258 mystream.close() 259 260 # Print to stdout 261 calc.txt = '-' 262 with calc.txt_outstream() as out: 263 calc.run(out=out) # output is written to stdout 264 """ 265 266 txt = self.txt 267 open_and_close = False # Do we open the file? 268 269 if txt is None: 270 # Suppress stdout 271 out = subprocess.DEVNULL 272 else: 273 if isinstance(txt, str): 274 if txt == '-': 275 # subprocess.call redirects this to stdout 276 out = None 277 else: 278 # Open the file in the work directory 279 txt = self._indir(txt) 280 # We wait with opening the file, until we are inside the 281 # try/finally 282 open_and_close = True 283 elif hasattr(txt, 'write'): 284 out = txt 285 else: 286 raise RuntimeError('txt should either be a string' 287 'or an I/O stream, got {}'.format(txt)) 288 289 try: 290 if open_and_close: 291 out = open(txt, 'w') 292 yield out 293 finally: 294 if open_and_close: 295 out.close() 296 297 def calculate(self, 298 atoms=None, 299 properties=('energy', ), 300 system_changes=tuple(calculator.all_changes)): 301 """Do a VASP calculation in the specified directory. 302 303 This will generate the necessary VASP input files, and then 304 execute VASP. After execution, the energy, forces. etc. are read 305 from the VASP output files. 306 """ 307 # Check for zero-length lattice vectors and PBC 308 # and that we actually have an Atoms object. 309 check_atoms(atoms) 310 311 self.clear_results() 312 313 if atoms is not None: 314 self.atoms = atoms.copy() 315 316 command = self.make_command(self.command) 317 self.write_input(self.atoms, properties, system_changes) 318 319 with self._txt_outstream() as out: 320 errorcode = self._run(command=command, 321 out=out, 322 directory=self.directory) 323 324 if errorcode: 325 raise calculator.CalculationFailed( 326 '{} in {} returned an error: {:d}'.format( 327 self.name, self.directory, errorcode)) 328 329 # Read results from calculation 330 self.update_atoms(atoms) 331 self.read_results() 332 333 def _run(self, command=None, out=None, directory=None): 334 """Method to explicitly execute VASP""" 335 if command is None: 336 command = self.command 337 if directory is None: 338 directory = self.directory 339 errorcode = subprocess.call(command, 340 shell=True, 341 stdout=out, 342 cwd=directory) 343 return errorcode 344 345 def check_state(self, atoms, tol=1e-15): 346 """Check for system changes since last calculation.""" 347 def compare_dict(d1, d2): 348 """Helper function to compare dictionaries""" 349 # Use symmetric difference to find keys which aren't shared 350 # for python 2.7 compatibility 351 if set(d1.keys()) ^ set(d2.keys()): 352 return False 353 354 # Check for differences in values 355 for key, value in d1.items(): 356 if np.any(value != d2[key]): 357 return False 358 return True 359 360 # First we check for default changes 361 system_changes = Calculator.check_state(self, atoms, tol=tol) 362 363 # We now check if we have made any changes to the input parameters 364 # XXX: Should we add these parameters to all_changes? 365 for param_string, old_dict in self.param_state.items(): 366 param_dict = getattr(self, param_string) # Get current param dict 367 if not compare_dict(param_dict, old_dict): 368 system_changes.append(param_string) 369 370 return system_changes 371 372 def _store_param_state(self): 373 """Store current parameter state""" 374 self.param_state = dict( 375 float_params=self.float_params.copy(), 376 exp_params=self.exp_params.copy(), 377 string_params=self.string_params.copy(), 378 int_params=self.int_params.copy(), 379 input_params=self.input_params.copy(), 380 bool_params=self.bool_params.copy(), 381 list_int_params=self.list_int_params.copy(), 382 list_bool_params=self.list_bool_params.copy(), 383 list_float_params=self.list_float_params.copy(), 384 dict_params=self.dict_params.copy()) 385 386 def asdict(self): 387 """Return a dictionary representation of the calculator state. 388 Does NOT contain information on the ``command``, ``txt`` or 389 ``directory`` keywords. 390 Contains the following keys: 391 392 - ``ase_version`` 393 - ``vasp_version`` 394 - ``inputs`` 395 - ``results`` 396 - ``atoms`` (Only if the calculator has an ``Atoms`` object) 397 """ 398 # Get versions 399 asevers = ase.__version__ 400 vaspvers = self.get_version() 401 402 self._store_param_state() # Update param state 403 # Store input parameters which have been set 404 inputs = { 405 key: value 406 for param_dct in self.param_state.values() 407 for key, value in param_dct.items() if value is not None 408 } 409 410 dct = { 411 'ase_version': asevers, 412 'vasp_version': vaspvers, 413 # '__ase_objtype__': self.ase_objtype, 414 'inputs': inputs, 415 'results': self.results.copy() 416 } 417 418 if self.atoms: 419 # Encode atoms as dict 420 from ase.db.row import atoms2dict 421 dct['atoms'] = atoms2dict(self.atoms) 422 423 return dct 424 425 def fromdict(self, dct): 426 """Restore calculator from a :func:`~ase.calculators.vasp.Vasp.asdict` 427 dictionary. 428 429 Parameters: 430 431 dct: Dictionary 432 The dictionary which is used to restore the calculator state. 433 """ 434 if 'vasp_version' in dct: 435 self.version = dct['vasp_version'] 436 if 'inputs' in dct: 437 self.set(**dct['inputs']) 438 self._store_param_state() 439 if 'atoms' in dct: 440 from ase.db.row import AtomsRow 441 atoms = AtomsRow(dct['atoms']).toatoms() 442 self.atoms = atoms 443 if 'results' in dct: 444 self.results.update(dct['results']) 445 446 def write_json(self, filename): 447 """Dump calculator state to JSON file. 448 449 Parameters: 450 451 filename: string 452 The filename which the JSON file will be stored to. 453 Prepends the ``directory`` path to the filename. 454 """ 455 filename = self._indir(filename) 456 dct = self.asdict() 457 jsonio.write_json(filename, dct) 458 459 def read_json(self, filename): 460 """Load Calculator state from an exported JSON Vasp file.""" 461 dct = jsonio.read_json(filename) 462 self.fromdict(dct) 463 464 def write_input(self, atoms, properties=None, system_changes=None): 465 """Write VASP inputfiles, INCAR, KPOINTS and POTCAR""" 466 # Create the folders where we write the files, if we aren't in the 467 # current working directory. 468 if self.directory != os.curdir and not os.path.isdir(self.directory): 469 os.makedirs(self.directory) 470 471 self.initialize(atoms) 472 473 GenerateVaspInput.write_input(self, atoms, directory=self.directory) 474 475 def read(self, label=None): 476 """Read results from VASP output files. 477 Files which are read: OUTCAR, CONTCAR and vasprun.xml 478 Raises ReadError if they are not found""" 479 if label is None: 480 label = self.label 481 Calculator.read(self, label) 482 483 # If we restart, self.parameters isn't initialized 484 if self.parameters is None: 485 self.parameters = self.get_default_parameters() 486 487 # Check for existence of the necessary output files 488 for f in ['OUTCAR', 'CONTCAR', 'vasprun.xml']: 489 file = self._indir(f) 490 if not file.is_file(): 491 raise calculator.ReadError( 492 'VASP outputfile {} was not found'.format(file)) 493 494 # Build sorting and resorting lists 495 self.read_sort() 496 497 # Read atoms 498 self.atoms = self.read_atoms(filename=self._indir('CONTCAR')) 499 500 # Read parameters 501 self.read_incar(filename=self._indir('INCAR')) 502 self.read_kpoints(filename=self._indir('KPOINTS')) 503 self.read_potcar(filename=self._indir('POTCAR')) 504 505 # Read the results from the calculation 506 self.read_results() 507 508 def _indir(self, filename): 509 """Prepend current directory to filename""" 510 return Path(self.directory) / filename 511 512 def read_sort(self): 513 """Create the sorting and resorting list from ase-sort.dat. 514 If the ase-sort.dat file does not exist, the sorting is redone. 515 """ 516 sortfile = self._indir('ase-sort.dat') 517 if os.path.isfile(sortfile): 518 self.sort = [] 519 self.resort = [] 520 with open(sortfile, 'r') as fd: 521 for line in fd: 522 sort, resort = line.split() 523 self.sort.append(int(sort)) 524 self.resort.append(int(resort)) 525 else: 526 # Redo the sorting 527 atoms = read(self._indir('CONTCAR')) 528 self.initialize(atoms) 529 530 def read_atoms(self, filename): 531 """Read the atoms from file located in the VASP 532 working directory. Normally called CONTCAR.""" 533 return read(filename)[self.resort] 534 535 def update_atoms(self, atoms): 536 """Update the atoms object with new positions and cell""" 537 if (self.int_params['ibrion'] is not None 538 and self.int_params['nsw'] is not None): 539 if self.int_params['ibrion'] > -1 and self.int_params['nsw'] > 0: 540 # Update atomic positions and unit cell with the ones read 541 # from CONTCAR. 542 atoms_sorted = read(self._indir('CONTCAR')) 543 atoms.positions = atoms_sorted[self.resort].positions 544 atoms.cell = atoms_sorted.cell 545 546 self.atoms = atoms # Creates a copy 547 548 def read_results(self): 549 """Read the results from VASP output files""" 550 # Temporarily load OUTCAR into memory 551 outcar = self.load_file('OUTCAR') 552 553 # Read the data we can from vasprun.xml 554 calc_xml = self._read_xml() 555 xml_results = calc_xml.results 556 557 # Fix sorting 558 xml_results['forces'] = xml_results['forces'][self.resort] 559 560 self.results.update(xml_results) 561 562 # Parse the outcar, as some properties are not loaded in vasprun.xml 563 # We want to limit this as much as possible, as reading large OUTCAR's 564 # is relatively slow 565 # Removed for now 566 # self.read_outcar(lines=outcar) 567 568 # Update results dict with results from OUTCAR 569 # which aren't written to the atoms object we read from 570 # the vasprun.xml file. 571 572 self.converged = self.read_convergence(lines=outcar) 573 self.version = self.read_version() 574 magmom, magmoms = self.read_mag(lines=outcar) 575 dipole = self.read_dipole(lines=outcar) 576 nbands = self.read_nbands(lines=outcar) 577 self.results.update( 578 dict(magmom=magmom, magmoms=magmoms, dipole=dipole, nbands=nbands)) 579 580 # Stress is not always present. 581 # Prevent calculation from going into a loop 582 if 'stress' not in self.results: 583 self.results.update(dict(stress=None)) 584 585 self._set_old_keywords() 586 587 # Store the parameters used for this calculation 588 self._store_param_state() 589 590 def _set_old_keywords(self): 591 """Store keywords for backwards compatibility wd VASP calculator""" 592 self.spinpol = self.get_spin_polarized() 593 self.energy_free = self.get_potential_energy(force_consistent=True) 594 self.energy_zero = self.get_potential_energy(force_consistent=False) 595 self.forces = self.get_forces() 596 self.fermi = self.get_fermi_level() 597 self.dipole = self.get_dipole_moment() 598 # Prevent calculation from going into a loop 599 self.stress = self.get_property('stress', allow_calculation=False) 600 self.nbands = self.get_number_of_bands() 601 602 # Below defines some functions for faster access to certain common keywords 603 @property 604 def kpts(self): 605 """Access the kpts from input_params dict""" 606 return self.input_params['kpts'] 607 608 @kpts.setter 609 def kpts(self, kpts): 610 """Set kpts in input_params dict""" 611 self.input_params['kpts'] = kpts 612 613 @property 614 def encut(self): 615 """Direct access to the encut parameter""" 616 return self.float_params['encut'] 617 618 @encut.setter 619 def encut(self, encut): 620 """Direct access for setting the encut parameter""" 621 self.set(encut=encut) 622 623 @property 624 def xc(self): 625 """Direct access to the xc parameter""" 626 return self.get_xc_functional() 627 628 @xc.setter 629 def xc(self, xc): 630 """Direct access for setting the xc parameter""" 631 self.set(xc=xc) 632 633 @property 634 def atoms(self): 635 return self._atoms 636 637 @atoms.setter 638 def atoms(self, atoms): 639 if atoms is None: 640 self._atoms = None 641 self.clear_results() 642 else: 643 if self.check_state(atoms): 644 self.clear_results() 645 self._atoms = atoms.copy() 646 647 def load_file(self, filename): 648 """Reads a file in the directory, and returns the lines 649 650 Example: 651 >>> outcar = load_file('OUTCAR') 652 """ 653 filename = self._indir(filename) 654 with open(filename, 'r') as fd: 655 return fd.readlines() 656 657 @contextmanager 658 def load_file_iter(self, filename): 659 """Return a file iterator""" 660 661 filename = self._indir(filename) 662 with open(filename, 'r') as fd: 663 yield fd 664 665 def read_outcar(self, lines=None): 666 """Read results from the OUTCAR file. 667 Deprecated, see read_results()""" 668 if not lines: 669 lines = self.load_file('OUTCAR') 670 # Spin polarized calculation? 671 self.spinpol = self.get_spin_polarized() 672 673 self.version = self.get_version() 674 675 # XXX: Do we want to read all of this again? 676 self.energy_free, self.energy_zero = self.read_energy(lines=lines) 677 self.forces = self.read_forces(lines=lines) 678 self.fermi = self.read_fermi(lines=lines) 679 680 self.dipole = self.read_dipole(lines=lines) 681 682 self.stress = self.read_stress(lines=lines) 683 self.nbands = self.read_nbands(lines=lines) 684 685 self.read_ldau() 686 self.magnetic_moment, self.magnetic_moments = self.read_mag( 687 lines=lines) 688 689 def _read_xml(self) -> SinglePointDFTCalculator: 690 """Read vasprun.xml, and return the last calculator object. 691 Returns calculator from the xml file. 692 Raises a ReadError if the reader is not able to construct a calculator. 693 """ 694 file = self._indir('vasprun.xml') 695 incomplete_msg = ( 696 f'The file "{file}" is incomplete, and no DFT data was available. ' 697 'This is likely due to an incomplete calculation.') 698 try: 699 _xml_atoms = read(file, index=-1, format='vasp-xml') 700 # Silence mypy, we should only ever get a single atoms object 701 assert isinstance(_xml_atoms, ase.Atoms) 702 except ElementTree.ParseError as exc: 703 raise calculator.ReadError(incomplete_msg) from exc 704 705 if _xml_atoms is None or _xml_atoms.calc is None: 706 raise calculator.ReadError(incomplete_msg) 707 708 self._xml_calc = _xml_atoms.calc 709 return self._xml_calc 710 711 @property 712 def _xml_calc(self) -> SinglePointDFTCalculator: 713 if self.__xml_calc is None: 714 raise RuntimeError(('vasprun.xml data has not yet been loaded. ' 715 'Run read_results() first.')) 716 return self.__xml_calc 717 718 @_xml_calc.setter 719 def _xml_calc(self, value): 720 self.__xml_calc = value 721 722 def get_ibz_k_points(self): 723 calc = self._xml_calc 724 return calc.get_ibz_k_points() 725 726 def get_kpt(self, kpt=0, spin=0): 727 calc = self._xml_calc 728 return calc.get_kpt(kpt=kpt, spin=spin) 729 730 def get_eigenvalues(self, kpt=0, spin=0): 731 calc = self._xml_calc 732 return calc.get_eigenvalues(kpt=kpt, spin=spin) 733 734 def get_fermi_level(self): 735 calc = self._xml_calc 736 return calc.get_fermi_level() 737 738 def get_homo_lumo(self): 739 calc = self._xml_calc 740 return calc.get_homo_lumo() 741 742 def get_homo_lumo_by_spin(self, spin=0): 743 calc = self._xml_calc 744 return calc.get_homo_lumo_by_spin(spin=spin) 745 746 def get_occupation_numbers(self, kpt=0, spin=0): 747 calc = self._xml_calc 748 return calc.get_occupation_numbers(kpt, spin) 749 750 def get_spin_polarized(self): 751 calc = self._xml_calc 752 return calc.get_spin_polarized() 753 754 def get_number_of_spins(self): 755 calc = self._xml_calc 756 return calc.get_number_of_spins() 757 758 def get_number_of_bands(self): 759 return self.results.get('nbands', None) 760 761 def get_number_of_electrons(self, lines=None): 762 if not lines: 763 lines = self.load_file('OUTCAR') 764 765 nelect = None 766 for line in lines: 767 if 'total number of electrons' in line: 768 nelect = float(line.split('=')[1].split()[0].strip()) 769 break 770 return nelect 771 772 def get_k_point_weights(self): 773 filename = self._indir('IBZKPT') 774 return self.read_k_point_weights(filename) 775 776 def get_dos(self, spin=None, **kwargs): 777 """ 778 The total DOS. 779 780 Uses the ASE DOS module, and returns a tuple with 781 (energies, dos). 782 """ 783 from ase.dft.dos import DOS 784 dos = DOS(self, **kwargs) 785 e = dos.get_energies() 786 d = dos.get_dos(spin=spin) 787 return e, d 788 789 def get_version(self): 790 if self.version is None: 791 # Try if we can read the version number 792 self.version = self.read_version() 793 return self.version 794 795 def read_version(self): 796 """Get the VASP version number""" 797 # The version number is the first occurrence, so we can just 798 # load the OUTCAR, as we will return soon anyway 799 if not os.path.isfile(self._indir('OUTCAR')): 800 return None 801 with self.load_file_iter('OUTCAR') as lines: 802 for line in lines: 803 if ' vasp.' in line: 804 return line[len(' vasp.'):].split()[0] 805 # We didn't find the version in VASP 806 return None 807 808 def get_number_of_iterations(self): 809 return self.read_number_of_iterations() 810 811 def read_number_of_iterations(self): 812 niter = None 813 with self.load_file_iter('OUTCAR') as lines: 814 for line in lines: 815 # find the last iteration number 816 if '- Iteration' in line: 817 niter = list(map(int, re.findall(r'\d+', line)))[1] 818 return niter 819 820 def read_number_of_ionic_steps(self): 821 niter = None 822 with self.load_file_iter('OUTCAR') as lines: 823 for line in lines: 824 if '- Iteration' in line: 825 niter = list(map(int, re.findall(r'\d+', line)))[0] 826 return niter 827 828 def read_stress(self, lines=None): 829 """Read stress from OUTCAR. 830 831 Depreciated: Use get_stress() instead. 832 """ 833 # We don't really need this, as we read this from vasprun.xml 834 # keeping it around "just in case" for now 835 if not lines: 836 lines = self.load_file('OUTCAR') 837 838 stress = None 839 for line in lines: 840 if ' in kB ' in line: 841 stress = -np.array([float(a) for a in line.split()[2:]]) 842 stress = stress[[0, 1, 2, 4, 5, 3]] * 1e-1 * ase.units.GPa 843 return stress 844 845 def read_ldau(self, lines=None): 846 """Read the LDA+U values from OUTCAR""" 847 if not lines: 848 lines = self.load_file('OUTCAR') 849 850 ldau_luj = None 851 ldauprint = None 852 ldau = None 853 ldautype = None 854 atomtypes = [] 855 # read ldau parameters from outcar 856 for line in lines: 857 if line.find('TITEL') != -1: # What atoms are present 858 atomtypes.append(line.split()[3].split('_')[0].split('.')[0]) 859 if line.find('LDAUTYPE') != -1: # Is this a DFT+U calculation 860 ldautype = int(line.split('=')[-1]) 861 ldau = True 862 ldau_luj = {} 863 if line.find('LDAUL') != -1: 864 L = line.split('=')[-1].split() 865 if line.find('LDAUU') != -1: 866 U = line.split('=')[-1].split() 867 if line.find('LDAUJ') != -1: 868 J = line.split('=')[-1].split() 869 # create dictionary 870 if ldau: 871 for i, symbol in enumerate(atomtypes): 872 ldau_luj[symbol] = { 873 'L': int(L[i]), 874 'U': float(U[i]), 875 'J': float(J[i]) 876 } 877 self.dict_params['ldau_luj'] = ldau_luj 878 879 self.ldau = ldau 880 self.ldauprint = ldauprint 881 self.ldautype = ldautype 882 self.ldau_luj = ldau_luj 883 return ldau, ldauprint, ldautype, ldau_luj 884 885 def get_xc_functional(self): 886 """Returns the XC functional or the pseudopotential type 887 888 If a XC recipe is set explicitly with 'xc', this is returned. 889 Otherwise, the XC functional associated with the 890 pseudopotentials (LDA, PW91 or PBE) is returned. 891 The string is always cast to uppercase for consistency 892 in checks.""" 893 if self.input_params.get('xc', None): 894 return self.input_params['xc'].upper() 895 if self.input_params.get('pp', None): 896 return self.input_params['pp'].upper() 897 raise ValueError('No xc or pp found.') 898 899 # Methods for reading information from OUTCAR files: 900 def read_energy(self, all=None, lines=None): 901 """Method to read energy from OUTCAR file. 902 Depreciated: use get_potential_energy() instead""" 903 if not lines: 904 lines = self.load_file('OUTCAR') 905 906 [energy_free, energy_zero] = [0, 0] 907 if all: 908 energy_free = [] 909 energy_zero = [] 910 for line in lines: 911 # Free energy 912 if line.lower().startswith(' free energy toten'): 913 if all: 914 energy_free.append(float(line.split()[-2])) 915 else: 916 energy_free = float(line.split()[-2]) 917 # Extrapolated zero point energy 918 if line.startswith(' energy without entropy'): 919 if all: 920 energy_zero.append(float(line.split()[-1])) 921 else: 922 energy_zero = float(line.split()[-1]) 923 return [energy_free, energy_zero] 924 925 def read_forces(self, all=False, lines=None): 926 """Method that reads forces from OUTCAR file. 927 928 If 'all' is switched on, the forces for all ionic steps 929 in the OUTCAR file be returned, in other case only the 930 forces for the last ionic configuration is returned.""" 931 932 if not lines: 933 lines = self.load_file('OUTCAR') 934 935 if all: 936 all_forces = [] 937 938 for n, line in enumerate(lines): 939 if 'TOTAL-FORCE' in line: 940 forces = [] 941 for i in range(len(self.atoms)): 942 forces.append( 943 np.array( 944 [float(f) for f in lines[n + 2 + i].split()[3:6]])) 945 946 if all: 947 all_forces.append(np.array(forces)[self.resort]) 948 949 if all: 950 return np.array(all_forces) 951 return np.array(forces)[self.resort] 952 953 def read_fermi(self, lines=None): 954 """Method that reads Fermi energy from OUTCAR file""" 955 if not lines: 956 lines = self.load_file('OUTCAR') 957 958 E_f = None 959 for line in lines: 960 if 'E-fermi' in line: 961 E_f = float(line.split()[2]) 962 return E_f 963 964 def read_dipole(self, lines=None): 965 """Read dipole from OUTCAR""" 966 if not lines: 967 lines = self.load_file('OUTCAR') 968 969 dipolemoment = np.zeros([1, 3]) 970 for line in lines: 971 if 'dipolmoment' in line: 972 dipolemoment = np.array([float(f) for f in line.split()[1:4]]) 973 return dipolemoment 974 975 def read_mag(self, lines=None): 976 if not lines: 977 lines = self.load_file('OUTCAR') 978 p = self.int_params 979 q = self.list_float_params 980 if self.spinpol: 981 magnetic_moment = self._read_magnetic_moment(lines=lines) 982 if ((p['lorbit'] is not None and p['lorbit'] >= 10) 983 or (p['lorbit'] is None and q['rwigs'])): 984 magnetic_moments = self._read_magnetic_moments(lines=lines) 985 else: 986 warn(('Magnetic moment data not written in OUTCAR (LORBIT<10),' 987 ' setting magnetic_moments to zero.\nSet LORBIT>=10' 988 ' to get information on magnetic moments')) 989 magnetic_moments = np.zeros(len(self.atoms)) 990 else: 991 magnetic_moment = 0.0 992 magnetic_moments = np.zeros(len(self.atoms)) 993 return magnetic_moment, magnetic_moments 994 995 def _read_magnetic_moments(self, lines=None): 996 """Read magnetic moments from OUTCAR. 997 Only reads the last occurrence. """ 998 if not lines: 999 lines = self.load_file('OUTCAR') 1000 1001 magnetic_moments = np.zeros(len(self.atoms)) 1002 magstr = 'magnetization (x)' 1003 1004 # Search for the last occurrence 1005 nidx = -1 1006 for n, line in enumerate(lines): 1007 if magstr in line: 1008 nidx = n 1009 1010 # Read that occurrence 1011 if nidx > -1: 1012 for m in range(len(self.atoms)): 1013 magnetic_moments[m] = float(lines[nidx + m + 4].split()[4]) 1014 return magnetic_moments[self.resort] 1015 1016 def _read_magnetic_moment(self, lines=None): 1017 """Read magnetic moment from OUTCAR""" 1018 if not lines: 1019 lines = self.load_file('OUTCAR') 1020 1021 for n, line in enumerate(lines): 1022 if 'number of electron ' in line: 1023 magnetic_moment = float(line.split()[-1]) 1024 return magnetic_moment 1025 1026 def read_nbands(self, lines=None): 1027 """Read number of bands from OUTCAR""" 1028 if not lines: 1029 lines = self.load_file('OUTCAR') 1030 1031 for line in lines: 1032 line = self.strip_warnings(line) 1033 if 'NBANDS' in line: 1034 return int(line.split()[-1]) 1035 return None 1036 1037 def read_convergence(self, lines=None): 1038 """Method that checks whether a calculation has converged.""" 1039 if not lines: 1040 lines = self.load_file('OUTCAR') 1041 1042 converged = None 1043 # First check electronic convergence 1044 for line in lines: 1045 if 0: # vasp always prints that! 1046 if line.rfind('aborting loop') > -1: # scf failed 1047 raise RuntimeError(line.strip()) 1048 break 1049 if 'EDIFF ' in line: 1050 ediff = float(line.split()[2]) 1051 if 'total energy-change' in line: 1052 # I saw this in an atomic oxygen calculation. it 1053 # breaks this code, so I am checking for it here. 1054 if 'MIXING' in line: 1055 continue 1056 split = line.split(':') 1057 a = float(split[1].split('(')[0]) 1058 b = split[1].split('(')[1][0:-2] 1059 # sometimes this line looks like (second number wrong format!): 1060 # energy-change (2. order) :-0.2141803E-08 ( 0.2737684-111) 1061 # we are checking still the first number so 1062 # let's "fix" the format for the second one 1063 if 'e' not in b.lower(): 1064 # replace last occurrence of - (assumed exponent) with -e 1065 bsplit = b.split('-') 1066 bsplit[-1] = 'e' + bsplit[-1] 1067 b = '-'.join(bsplit).replace('-e', 'e-') 1068 b = float(b) 1069 if [abs(a), abs(b)] < [ediff, ediff]: 1070 converged = True 1071 else: 1072 converged = False 1073 continue 1074 # Then if ibrion in [1,2,3] check whether ionic relaxation 1075 # condition been fulfilled 1076 if ((self.int_params['ibrion'] in [1, 2, 3] 1077 and self.int_params['nsw'] not in [0])): 1078 if not self.read_relaxed(): 1079 converged = False 1080 else: 1081 converged = True 1082 return converged 1083 1084 def read_k_point_weights(self, filename): 1085 """Read k-point weighting. Normally named IBZKPT.""" 1086 1087 lines = self.load_file(filename) 1088 1089 if 'Tetrahedra\n' in lines: 1090 N = lines.index('Tetrahedra\n') 1091 else: 1092 N = len(lines) 1093 kpt_weights = [] 1094 for n in range(3, N): 1095 kpt_weights.append(float(lines[n].split()[3])) 1096 kpt_weights = np.array(kpt_weights) 1097 kpt_weights /= np.sum(kpt_weights) 1098 1099 return kpt_weights 1100 1101 def read_relaxed(self, lines=None): 1102 """Check if ionic relaxation completed""" 1103 if not lines: 1104 lines = self.load_file('OUTCAR') 1105 for line in lines: 1106 if 'reached required accuracy' in line: 1107 return True 1108 return False 1109 1110 def read_spinpol(self, lines=None): 1111 """Method which reads if a calculation from spinpolarized using OUTCAR. 1112 1113 Depreciated: Use get_spin_polarized() instead. 1114 """ 1115 if not lines: 1116 lines = self.load_file('OUTCAR') 1117 1118 for line in lines: 1119 if 'ISPIN' in line: 1120 if int(line.split()[2]) == 2: 1121 self.spinpol = True 1122 else: 1123 self.spinpol = False 1124 return self.spinpol 1125 1126 def strip_warnings(self, line): 1127 """Returns empty string instead of line from warnings in OUTCAR.""" 1128 if line[0] == "|": 1129 return "" 1130 return line 1131 1132 @property 1133 def txt(self): 1134 return self._txt 1135 1136 @txt.setter 1137 def txt(self, txt): 1138 if isinstance(txt, PurePath): 1139 txt = str(txt) 1140 self._txt = txt 1141 1142 def get_number_of_grid_points(self): 1143 raise NotImplementedError 1144 1145 def get_pseudo_density(self): 1146 raise NotImplementedError 1147 1148 def get_pseudo_wavefunction(self, n=0, k=0, s=0, pad=True): 1149 raise NotImplementedError 1150 1151 def get_bz_k_points(self): 1152 raise NotImplementedError 1153 1154 def read_vib_freq(self, lines=None): 1155 """Read vibrational frequencies. 1156 1157 Returns list of real and list of imaginary frequencies.""" 1158 freq = [] 1159 i_freq = [] 1160 1161 if not lines: 1162 lines = self.load_file('OUTCAR') 1163 1164 for line in lines: 1165 data = line.split() 1166 if 'THz' in data: 1167 if 'f/i=' not in data: 1168 freq.append(float(data[-2])) 1169 else: 1170 i_freq.append(float(data[-2])) 1171 return freq, i_freq 1172 1173 def get_nonselfconsistent_energies(self, bee_type): 1174 """ Method that reads and returns BEE energy contributions 1175 written in OUTCAR file. 1176 """ 1177 assert bee_type == 'beefvdw' 1178 cmd = 'grep -32 "BEEF xc energy contributions" OUTCAR | tail -32' 1179 p = os.popen(cmd, 'r') 1180 s = p.readlines() 1181 p.close() 1182 xc = np.array([]) 1183 for line in s: 1184 l_ = float(line.split(":")[-1]) 1185 xc = np.append(xc, l_) 1186 assert len(xc) == 32 1187 return xc 1188 1189 1190##################################### 1191# Below defines helper functions 1192# for the VASP calculator 1193##################################### 1194 1195 1196def check_atoms(atoms: ase.Atoms) -> None: 1197 """Perform checks on the atoms object, to verify that 1198 it can be run by VASP. 1199 A CalculatorSetupError error is raised if the atoms are not supported. 1200 """ 1201 1202 # Loop through all check functions 1203 for check in (check_atoms_type, check_cell, check_pbc): 1204 check(atoms) 1205 1206 1207def check_cell(atoms: ase.Atoms) -> None: 1208 """Check if there is a zero unit cell. 1209 Raises CalculatorSetupError if the cell is wrong. 1210 """ 1211 if atoms.cell.rank < 3: 1212 raise calculator.CalculatorSetupError( 1213 "The lattice vectors are zero! " 1214 "This is the default value - please specify a " 1215 "unit cell.") 1216 1217 1218def check_pbc(atoms: ase.Atoms) -> None: 1219 """Check if any boundaries are not PBC, as VASP 1220 cannot handle non-PBC. 1221 Raises CalculatorSetupError. 1222 """ 1223 if not atoms.pbc.all(): 1224 raise calculator.CalculatorSetupError( 1225 "Vasp cannot handle non-periodic boundaries. " 1226 "Please enable all PBC, e.g. atoms.pbc=True") 1227 1228 1229def check_atoms_type(atoms: ase.Atoms) -> None: 1230 """Check that the passed atoms object is in fact an Atoms object. 1231 Raises CalculatorSetupError. 1232 """ 1233 if not isinstance(atoms, ase.Atoms): 1234 raise calculator.CalculatorSetupError( 1235 ('Expected an Atoms object, ' 1236 'instead got object of type {}'.format(type(atoms)))) 1237