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