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