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