Coverage for /builds/debichem-team/python-ase/ase/calculators/vasp/vasp.py: 41.42%
676 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-03-06 04:00 +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 Calculator.calculate(self, atoms, properties, system_changes)
325 # Check for zero-length lattice vectors and PBC
326 # and that we actually have an Atoms object.
327 check_atoms(self.atoms)
329 self.clear_results()
331 command = self.make_command(self.command)
332 self.write_input(self.atoms, properties, system_changes)
334 with self._txt_outstream() as out:
335 errorcode, stderr = self._run(command=command,
336 out=out,
337 directory=self.directory)
339 if errorcode:
340 raise calculator.CalculationFailed(
341 '{} in {} returned an error: {:d} stderr {}'.format(
342 self.name, Path(self.directory).resolve(), errorcode,
343 stderr))
345 # Read results from calculation
346 self.update_atoms(atoms)
347 self.read_results()
349 def _run(self, command=None, out=None, directory=None):
350 """Method to explicitly execute VASP"""
351 if command is None:
352 command = self.command
353 if directory is None:
354 directory = self.directory
355 result = subprocess.run(command,
356 shell=True,
357 cwd=directory,
358 capture_output=True,
359 text=True)
360 if out is not None:
361 out.write(result.stdout)
362 return result.returncode, result.stderr
364 def check_state(self, atoms, tol=1e-15):
365 """Check for system changes since last calculation."""
366 def compare_dict(d1, d2):
367 """Helper function to compare dictionaries"""
368 # Use symmetric difference to find keys which aren't shared
369 # for python 2.7 compatibility
370 if set(d1.keys()) ^ set(d2.keys()):
371 return False
373 # Check for differences in values
374 for key, value in d1.items():
375 if np.any(value != d2[key]):
376 return False
377 return True
379 # First we check for default changes
380 system_changes = Calculator.check_state(self, atoms, tol=tol)
382 # We now check if we have made any changes to the input parameters
383 # XXX: Should we add these parameters to all_changes?
384 for param_string, old_dict in self.param_state.items():
385 param_dict = getattr(self, param_string) # Get current param dict
386 if not compare_dict(param_dict, old_dict):
387 system_changes.append(param_string)
389 return system_changes
391 def _store_param_state(self):
392 """Store current parameter state"""
393 self.param_state = dict(
394 float_params=self.float_params.copy(),
395 exp_params=self.exp_params.copy(),
396 string_params=self.string_params.copy(),
397 int_params=self.int_params.copy(),
398 input_params=self.input_params.copy(),
399 bool_params=self.bool_params.copy(),
400 list_int_params=self.list_int_params.copy(),
401 list_bool_params=self.list_bool_params.copy(),
402 list_float_params=self.list_float_params.copy(),
403 dict_params=self.dict_params.copy(),
404 special_params=self.special_params.copy())
406 def asdict(self):
407 """Return a dictionary representation of the calculator state.
408 Does NOT contain information on the ``command``, ``txt`` or
409 ``directory`` keywords.
410 Contains the following keys:
412 - ``ase_version``
413 - ``vasp_version``
414 - ``inputs``
415 - ``results``
416 - ``atoms`` (Only if the calculator has an ``Atoms`` object)
417 """
418 # Get versions
419 asevers = ase.__version__
420 vaspvers = self.get_version()
422 self._store_param_state() # Update param state
423 # Store input parameters which have been set
424 inputs = {
425 key: value
426 for param_dct in self.param_state.values()
427 for key, value in param_dct.items() if value is not None
428 }
430 dct = {
431 'ase_version': asevers,
432 'vasp_version': vaspvers,
433 # '__ase_objtype__': self.ase_objtype,
434 'inputs': inputs,
435 'results': self.results.copy()
436 }
438 if self.atoms:
439 # Encode atoms as dict
440 from ase.db.row import atoms2dict
441 dct['atoms'] = atoms2dict(self.atoms)
443 return dct
445 def fromdict(self, dct):
446 """Restore calculator from a :func:`~ase.calculators.vasp.Vasp.asdict`
447 dictionary.
449 Parameters:
451 dct: Dictionary
452 The dictionary which is used to restore the calculator state.
453 """
454 if 'vasp_version' in dct:
455 self.version = dct['vasp_version']
456 if 'inputs' in dct:
457 self.set(**dct['inputs'])
458 self._store_param_state()
459 if 'atoms' in dct:
460 from ase.db.row import AtomsRow
461 atoms = AtomsRow(dct['atoms']).toatoms()
462 self.atoms = atoms
463 if 'results' in dct:
464 self.results.update(dct['results'])
466 def write_json(self, filename):
467 """Dump calculator state to JSON file.
469 Parameters:
471 filename: string
472 The filename which the JSON file will be stored to.
473 Prepends the ``directory`` path to the filename.
474 """
475 filename = self._indir(filename)
476 dct = self.asdict()
477 jsonio.write_json(filename, dct)
479 def read_json(self, filename):
480 """Load Calculator state from an exported JSON Vasp file."""
481 dct = jsonio.read_json(filename)
482 self.fromdict(dct)
484 def write_input(self, atoms, properties=None, system_changes=None):
485 """Write VASP inputfiles, INCAR, KPOINTS and POTCAR"""
486 self.initialize(atoms)
487 GenerateVaspInput.write_input(self, atoms, directory=self.directory)
489 def read(self, label=None):
490 """Read results from VASP output files.
491 Files which are read: OUTCAR, CONTCAR and vasprun.xml
492 Raises ReadError if they are not found"""
493 if label is None:
494 label = self.label
495 Calculator.read(self, label)
497 # If we restart, self.parameters isn't initialized
498 if self.parameters is None:
499 self.parameters = self.get_default_parameters()
501 # Check for existence of the necessary output files
502 for f in ['OUTCAR', 'CONTCAR', 'vasprun.xml']:
503 file = self._indir(f)
504 if not file.is_file():
505 raise calculator.ReadError(
506 f'VASP outputfile {file} was not found')
508 # Build sorting and resorting lists
509 self.read_sort()
511 # Read atoms
512 self.atoms = self.read_atoms(filename=self._indir('CONTCAR'))
514 # Read parameters
515 self.read_incar(filename=self._indir('INCAR'))
516 self.read_kpoints(filename=self._indir('KPOINTS'))
517 self.read_potcar(filename=self._indir('POTCAR'))
519 # Read the results from the calculation
520 self.read_results()
522 def _indir(self, filename):
523 """Prepend current directory to filename"""
524 return Path(self.directory) / filename
526 def read_sort(self):
527 """Create the sorting and resorting list from ase-sort.dat.
528 If the ase-sort.dat file does not exist, the sorting is redone.
529 """
530 sortfile = self._indir('ase-sort.dat')
531 if os.path.isfile(sortfile):
532 self.sort = []
533 self.resort = []
534 with open(sortfile) as fd:
535 for line in fd:
536 sort, resort = line.split()
537 self.sort.append(int(sort))
538 self.resort.append(int(resort))
539 else:
540 # Redo the sorting
541 atoms = read(self._indir('CONTCAR'))
542 self.initialize(atoms)
544 def read_atoms(self, filename):
545 """Read the atoms from file located in the VASP
546 working directory. Normally called CONTCAR."""
547 return read(filename)[self.resort]
549 def update_atoms(self, atoms):
550 """Update the atoms object with new positions and cell"""
551 if (self.int_params['ibrion'] is not None
552 and self.int_params['nsw'] is not None):
553 if self.int_params['ibrion'] > -1 and self.int_params['nsw'] > 0:
554 # Update atomic positions and unit cell with the ones read
555 # from CONTCAR.
556 atoms_sorted = read(self._indir('CONTCAR'))
557 atoms.positions = atoms_sorted[self.resort].positions
558 atoms.cell = atoms_sorted.cell
560 self.atoms = atoms # Creates a copy
562 def read_results(self):
563 """Read the results from VASP output files"""
564 # Temporarily load OUTCAR into memory
565 outcar = self.load_file('OUTCAR')
567 # Read the data we can from vasprun.xml
568 calc_xml = self._read_xml()
569 xml_results = calc_xml.results
571 # Fix sorting
572 xml_results['forces'] = xml_results['forces'][self.resort]
574 self.results.update(xml_results)
576 # Parse the outcar, as some properties are not loaded in vasprun.xml
577 # We want to limit this as much as possible, as reading large OUTCAR's
578 # is relatively slow
579 # Removed for now
580 # self.read_outcar(lines=outcar)
582 # Update results dict with results from OUTCAR
583 # which aren't written to the atoms object we read from
584 # the vasprun.xml file.
586 self.converged = self.read_convergence(lines=outcar)
587 self.version = self.read_version()
588 magmom, magmoms = self.read_mag(lines=outcar)
589 dipole = self.read_dipole(lines=outcar)
590 nbands = self.read_nbands(lines=outcar)
591 self.results.update(
592 dict(magmom=magmom, magmoms=magmoms, dipole=dipole, nbands=nbands))
594 # Stress is not always present.
595 # Prevent calculation from going into a loop
596 if 'stress' not in self.results:
597 self.results.update(dict(stress=None))
599 self._set_old_keywords()
601 # Store the parameters used for this calculation
602 self._store_param_state()
604 def _set_old_keywords(self):
605 """Store keywords for backwards compatibility wd VASP calculator"""
606 self.spinpol = self.get_spin_polarized()
607 self.energy_free = self.get_potential_energy(force_consistent=True)
608 self.energy_zero = self.get_potential_energy(force_consistent=False)
609 self.forces = self.get_forces()
610 self.fermi = self.get_fermi_level()
611 self.dipole = self.get_dipole_moment()
612 # Prevent calculation from going into a loop
613 self.stress = self.get_property('stress', allow_calculation=False)
614 self.nbands = self.get_number_of_bands()
616 # Below defines some functions for faster access to certain common keywords
617 @property
618 def kpts(self):
619 """Access the kpts from input_params dict"""
620 return self.input_params['kpts']
622 @kpts.setter
623 def kpts(self, kpts):
624 """Set kpts in input_params dict"""
625 self.input_params['kpts'] = kpts
627 @property
628 def encut(self):
629 """Direct access to the encut parameter"""
630 return self.float_params['encut']
632 @encut.setter
633 def encut(self, encut):
634 """Direct access for setting the encut parameter"""
635 self.set(encut=encut)
637 @property
638 def xc(self):
639 """Direct access to the xc parameter"""
640 return self.get_xc_functional()
642 @xc.setter
643 def xc(self, xc):
644 """Direct access for setting the xc parameter"""
645 self.set(xc=xc)
647 @property
648 def atoms(self):
649 return self._atoms
651 @atoms.setter
652 def atoms(self, atoms):
653 if atoms is None:
654 self._atoms = None
655 self.clear_results()
656 else:
657 if self.check_state(atoms):
658 self.clear_results()
659 self._atoms = atoms.copy()
661 def load_file(self, filename):
662 """Reads a file in the directory, and returns the lines
664 Example:
665 >>> from ase.calculators.vasp import Vasp
666 >>> calc = Vasp()
667 >>> outcar = calc.load_file('OUTCAR') # doctest: +SKIP
668 """
669 filename = self._indir(filename)
670 with open(filename) as fd:
671 return fd.readlines()
673 @contextmanager
674 def load_file_iter(self, filename):
675 """Return a file iterator"""
677 filename = self._indir(filename)
678 with open(filename) as fd:
679 yield fd
681 def read_outcar(self, lines=None):
682 """Read results from the OUTCAR file.
683 Deprecated, see read_results()"""
684 if not lines:
685 lines = self.load_file('OUTCAR')
686 # Spin polarized calculation?
687 self.spinpol = self.get_spin_polarized()
689 self.version = self.get_version()
691 # XXX: Do we want to read all of this again?
692 self.energy_free, self.energy_zero = self.read_energy(lines=lines)
693 self.forces = self.read_forces(lines=lines)
694 self.fermi = self.read_fermi(lines=lines)
696 self.dipole = self.read_dipole(lines=lines)
698 self.stress = self.read_stress(lines=lines)
699 self.nbands = self.read_nbands(lines=lines)
701 self.read_ldau()
702 self.magnetic_moment, self.magnetic_moments = self.read_mag(
703 lines=lines)
705 def _read_xml(self) -> SinglePointDFTCalculator:
706 """Read vasprun.xml, and return the last calculator object.
707 Returns calculator from the xml file.
708 Raises a ReadError if the reader is not able to construct a calculator.
709 """
710 file = self._indir('vasprun.xml')
711 incomplete_msg = (
712 f'The file "{file}" is incomplete, and no DFT data was available. '
713 'This is likely due to an incomplete calculation.')
714 try:
715 _xml_atoms = read(file, index=-1, format='vasp-xml')
716 # Silence mypy, we should only ever get a single atoms object
717 assert isinstance(_xml_atoms, ase.Atoms)
718 except ElementTree.ParseError as exc:
719 raise calculator.ReadError(incomplete_msg) from exc
721 if _xml_atoms is None or _xml_atoms.calc is None:
722 raise calculator.ReadError(incomplete_msg)
724 self._xml_calc = _xml_atoms.calc
725 return self._xml_calc
727 @property
728 def _xml_calc(self) -> SinglePointDFTCalculator:
729 if self.__xml_calc is None:
730 raise RuntimeError('vasprun.xml data has not yet been loaded. '
731 'Run read_results() first.')
732 return self.__xml_calc
734 @_xml_calc.setter
735 def _xml_calc(self, value):
736 self.__xml_calc = value
738 def get_ibz_k_points(self):
739 calc = self._xml_calc
740 return calc.get_ibz_k_points()
742 def get_kpt(self, kpt=0, spin=0):
743 calc = self._xml_calc
744 return calc.get_kpt(kpt=kpt, spin=spin)
746 def get_eigenvalues(self, kpt=0, spin=0):
747 calc = self._xml_calc
748 return calc.get_eigenvalues(kpt=kpt, spin=spin)
750 def get_fermi_level(self):
751 calc = self._xml_calc
752 return calc.get_fermi_level()
754 def get_homo_lumo(self):
755 calc = self._xml_calc
756 return calc.get_homo_lumo()
758 def get_homo_lumo_by_spin(self, spin=0):
759 calc = self._xml_calc
760 return calc.get_homo_lumo_by_spin(spin=spin)
762 def get_occupation_numbers(self, kpt=0, spin=0):
763 calc = self._xml_calc
764 return calc.get_occupation_numbers(kpt, spin)
766 def get_spin_polarized(self):
767 calc = self._xml_calc
768 return calc.get_spin_polarized()
770 def get_number_of_spins(self):
771 calc = self._xml_calc
772 return calc.get_number_of_spins()
774 def get_number_of_bands(self):
775 return self.results.get('nbands', None)
777 def get_number_of_electrons(self, lines=None):
778 if not lines:
779 lines = self.load_file('OUTCAR')
781 nelect = None
782 for line in lines:
783 if 'total number of electrons' in line:
784 nelect = float(line.split('=')[1].split()[0].strip())
785 break
786 return nelect
788 def get_k_point_weights(self):
789 filename = 'IBZKPT'
790 return self.read_k_point_weights(filename)
792 def get_dos(self, spin=None, **kwargs):
793 """
794 The total DOS.
796 Uses the ASE DOS module, and returns a tuple with
797 (energies, dos).
798 """
799 from ase.dft.dos import DOS
800 dos = DOS(self, **kwargs)
801 e = dos.get_energies()
802 d = dos.get_dos(spin=spin)
803 return e, d
805 def get_version(self):
806 if self.version is None:
807 # Try if we can read the version number
808 self.version = self.read_version()
809 return self.version
811 def read_version(self):
812 """Get the VASP version number"""
813 # The version number is the first occurrence, so we can just
814 # load the OUTCAR, as we will return soon anyway
815 if not os.path.isfile(self._indir('OUTCAR')):
816 return None
817 with self.load_file_iter('OUTCAR') as lines:
818 for line in lines:
819 if ' vasp.' in line:
820 return line[len(' vasp.'):].split()[0]
821 # We didn't find the version in VASP
822 return None
824 def get_number_of_iterations(self):
825 return self.read_number_of_iterations()
827 def read_number_of_iterations(self):
828 niter = None
829 with self.load_file_iter('OUTCAR') as lines:
830 for line in lines:
831 # find the last iteration number
832 if '- Iteration' in line:
833 niter = list(map(int, re.findall(r'\d+', line)))[1]
834 return niter
836 def read_number_of_ionic_steps(self):
837 niter = None
838 with self.load_file_iter('OUTCAR') as lines:
839 for line in lines:
840 if '- Iteration' in line:
841 niter = list(map(int, re.findall(r'\d+', line)))[0]
842 return niter
844 def read_stress(self, lines=None):
845 """Read stress from OUTCAR.
847 Depreciated: Use get_stress() instead.
848 """
849 # We don't really need this, as we read this from vasprun.xml
850 # keeping it around "just in case" for now
851 if not lines:
852 lines = self.load_file('OUTCAR')
854 stress = None
855 for line in lines:
856 if ' in kB ' in line:
857 stress = -np.array([float(a) for a in line.split()[2:]])
858 stress = stress[[0, 1, 2, 4, 5, 3]] * 1e-1 * ase.units.GPa
859 return stress
861 def read_ldau(self, lines=None):
862 """Read the LDA+U values from OUTCAR"""
863 if not lines:
864 lines = self.load_file('OUTCAR')
866 ldau_luj = None
867 ldauprint = None
868 ldau = None
869 ldautype = None
870 atomtypes = []
871 # read ldau parameters from outcar
872 for line in lines:
873 if line.find('TITEL') != -1: # What atoms are present
874 atomtypes.append(line.split()[3].split('_')[0].split('.')[0])
875 if line.find('LDAUTYPE') != -1: # Is this a DFT+U calculation
876 ldautype = int(line.split('=')[-1])
877 ldau = True
878 ldau_luj = {}
879 if line.find('LDAUL') != -1:
880 L = line.split('=')[-1].split()
881 if line.find('LDAUU') != -1:
882 U = line.split('=')[-1].split()
883 if line.find('LDAUJ') != -1:
884 J = line.split('=')[-1].split()
885 # create dictionary
886 if ldau:
887 for i, symbol in enumerate(atomtypes):
888 ldau_luj[symbol] = {
889 'L': int(L[i]),
890 'U': float(U[i]),
891 'J': float(J[i])
892 }
893 self.dict_params['ldau_luj'] = ldau_luj
895 self.ldau = ldau
896 self.ldauprint = ldauprint
897 self.ldautype = ldautype
898 self.ldau_luj = ldau_luj
899 return ldau, ldauprint, ldautype, ldau_luj
901 def get_xc_functional(self):
902 """Returns the XC functional or the pseudopotential type
904 If a XC recipe is set explicitly with 'xc', this is returned.
905 Otherwise, the XC functional associated with the
906 pseudopotentials (LDA, PW91 or PBE) is returned.
907 The string is always cast to uppercase for consistency
908 in checks."""
909 if self.input_params.get('xc', None):
910 return self.input_params['xc'].upper()
911 if self.input_params.get('pp', None):
912 return self.input_params['pp'].upper()
913 raise ValueError('No xc or pp found.')
915 # Methods for reading information from OUTCAR files:
916 def read_energy(self, all=None, lines=None):
917 """Method to read energy from OUTCAR file.
918 Depreciated: use get_potential_energy() instead"""
919 if not lines:
920 lines = self.load_file('OUTCAR')
922 [energy_free, energy_zero] = [0, 0]
923 if all:
924 energy_free = []
925 energy_zero = []
926 for line in lines:
927 # Free energy
928 if line.lower().startswith(' free energy toten'):
929 if all:
930 energy_free.append(float(line.split()[-2]))
931 else:
932 energy_free = float(line.split()[-2])
933 # Extrapolated zero point energy
934 if line.startswith(' energy without entropy'):
935 if all:
936 energy_zero.append(float(line.split()[-1]))
937 else:
938 energy_zero = float(line.split()[-1])
939 return [energy_free, energy_zero]
941 def read_forces(self, all=False, lines=None):
942 """Method that reads forces from OUTCAR file.
944 If 'all' is switched on, the forces for all ionic steps
945 in the OUTCAR file be returned, in other case only the
946 forces for the last ionic configuration is returned."""
948 if not lines:
949 lines = self.load_file('OUTCAR')
951 if all:
952 all_forces = []
954 for n, line in enumerate(lines):
955 if 'TOTAL-FORCE' in line:
956 forces = []
957 for i in range(len(self.atoms)):
958 forces.append(
959 np.array(
960 [float(f) for f in lines[n + 2 + i].split()[3:6]]))
962 if all:
963 all_forces.append(np.array(forces)[self.resort])
965 if all:
966 return np.array(all_forces)
967 return np.array(forces)[self.resort]
969 def read_fermi(self, lines=None):
970 """Method that reads Fermi energy from OUTCAR file"""
971 if not lines:
972 lines = self.load_file('OUTCAR')
974 E_f = None
975 for line in lines:
976 if 'E-fermi' in line:
977 E_f = float(line.split()[2])
978 return E_f
980 def read_dipole(self, lines=None):
981 """Read dipole from OUTCAR"""
982 if not lines:
983 lines = self.load_file('OUTCAR')
985 dipolemoment = np.zeros([1, 3])
986 for line in lines:
987 if 'dipolmoment' in line:
988 dipolemoment = np.array([float(f) for f in line.split()[1:4]])
989 return dipolemoment
991 def read_mag(self, lines=None):
992 if not lines:
993 lines = self.load_file('OUTCAR')
994 p = self.int_params
995 q = self.list_float_params
996 if self.spinpol:
997 magnetic_moment = self._read_magnetic_moment(lines=lines)
998 if ((p['lorbit'] is not None and p['lorbit'] >= 10)
999 or (p['lorbit'] is None and q['rwigs'])):
1000 magnetic_moments = self._read_magnetic_moments(lines=lines)
1001 else:
1002 warn('Magnetic moment data not written in OUTCAR (LORBIT<10),'
1003 ' setting magnetic_moments to zero.\nSet LORBIT>=10'
1004 ' to get information on magnetic moments')
1005 magnetic_moments = np.zeros(len(self.atoms))
1006 else:
1007 magnetic_moment = 0.0
1008 magnetic_moments = np.zeros(len(self.atoms))
1009 return magnetic_moment, magnetic_moments
1011 def _read_magnetic_moments(self, lines=None):
1012 """Read magnetic moments from OUTCAR.
1013 Only reads the last occurrence. """
1014 if not lines:
1015 lines = self.load_file('OUTCAR')
1017 magnetic_moments = np.zeros(len(self.atoms))
1018 magstr = 'magnetization (x)'
1020 # Search for the last occurrence
1021 nidx = -1
1022 for n, line in enumerate(lines):
1023 if magstr in line:
1024 nidx = n
1026 # Read that occurrence
1027 if nidx > -1:
1028 for m in range(len(self.atoms)):
1029 magnetic_moments[m] = float(lines[nidx + m + 4].split()[-1])
1030 return magnetic_moments[self.resort]
1032 def _read_magnetic_moment(self, lines=None):
1033 """Read magnetic moment from OUTCAR"""
1034 if not lines:
1035 lines = self.load_file('OUTCAR')
1037 for line in lines:
1038 if 'number of electron ' in line:
1039 _ = line.split()[-1]
1040 magnetic_moment = 0.0 if _ == "magnetization" else float(_)
1041 return magnetic_moment
1043 def read_nbands(self, lines=None):
1044 """Read number of bands from OUTCAR"""
1045 if not lines:
1046 lines = self.load_file('OUTCAR')
1048 for line in lines:
1049 line = self.strip_warnings(line)
1050 if 'NBANDS' in line:
1051 return int(line.split()[-1])
1052 return None
1054 def read_convergence(self, lines=None):
1055 """Method that checks whether a calculation has converged."""
1056 if not lines:
1057 lines = self.load_file('OUTCAR')
1059 converged = None
1060 # First check electronic convergence
1061 for line in lines:
1062 # VASP 6 actually labels loop exit reason
1063 if 'aborting loop' in line:
1064 converged = 'because EDIFF is reached' in line
1065 # NOTE: 'EDIFF was not reached (unconverged)'
1066 # and
1067 # 'because hard stop was set'
1068 # will result in unconverged
1069 break
1070 # determine convergence by attempting to reproduce VASP's
1071 # internal logic
1072 if 'EDIFF ' in line:
1073 ediff = float(line.split()[2])
1074 if 'total energy-change' in line:
1075 # I saw this in an atomic oxygen calculation. it
1076 # breaks this code, so I am checking for it here.
1077 if 'MIXING' in line:
1078 continue
1079 split = line.split(':')
1080 a = float(split[1].split('(')[0])
1081 b = split[1].split('(')[1][0:-2]
1082 # sometimes this line looks like (second number wrong format!):
1083 # energy-change (2. order) :-0.2141803E-08 ( 0.2737684-111)
1084 # we are checking still the first number so
1085 # let's "fix" the format for the second one
1086 if 'e' not in b.lower():
1087 # replace last occurrence of - (assumed exponent) with -e
1088 bsplit = b.split('-')
1089 bsplit[-1] = 'e' + bsplit[-1]
1090 b = '-'.join(bsplit).replace('-e', 'e-')
1091 b = float(b)
1092 if [abs(a), abs(b)] < [ediff, ediff]:
1093 converged = True
1094 else:
1095 converged = False
1096 continue
1097 # Then if ibrion in [1,2,3] check whether ionic relaxation
1098 # condition been fulfilled
1099 if (self.int_params['ibrion'] in [1, 2, 3]
1100 and self.int_params['nsw'] not in [0]):
1101 if not self.read_relaxed():
1102 converged = False
1103 else:
1104 converged = True
1105 return converged
1107 def read_k_point_weights(self, filename):
1108 """Read k-point weighting. Normally named IBZKPT."""
1110 lines = self.load_file(filename)
1112 if 'Tetrahedra\n' in lines:
1113 N = lines.index('Tetrahedra\n')
1114 else:
1115 N = len(lines)
1116 kpt_weights = []
1117 for n in range(3, N):
1118 kpt_weights.append(float(lines[n].split()[3]))
1119 kpt_weights = np.array(kpt_weights)
1120 kpt_weights /= np.sum(kpt_weights)
1122 return kpt_weights
1124 def read_relaxed(self, lines=None):
1125 """Check if ionic relaxation completed"""
1126 if not lines:
1127 lines = self.load_file('OUTCAR')
1128 for line in lines:
1129 if 'reached required accuracy' in line:
1130 return True
1131 return False
1133 def read_spinpol(self, lines=None):
1134 """Method which reads if a calculation from spinpolarized using OUTCAR.
1136 Depreciated: Use get_spin_polarized() instead.
1137 """
1138 if not lines:
1139 lines = self.load_file('OUTCAR')
1141 for line in lines:
1142 if 'ISPIN' in line:
1143 if int(line.split()[2]) == 2:
1144 self.spinpol = True
1145 else:
1146 self.spinpol = False
1147 return self.spinpol
1149 def strip_warnings(self, line):
1150 """Returns empty string instead of line from warnings in OUTCAR."""
1151 if line[0] == "|":
1152 return ""
1153 return line
1155 @property
1156 def txt(self):
1157 return self._txt
1159 @txt.setter
1160 def txt(self, txt):
1161 if isinstance(txt, PurePath):
1162 txt = str(txt)
1163 self._txt = txt
1165 def get_number_of_grid_points(self):
1166 raise NotImplementedError
1168 def get_pseudo_density(self):
1169 raise NotImplementedError
1171 def get_pseudo_wavefunction(self, n=0, k=0, s=0, pad=True):
1172 raise NotImplementedError
1174 def get_bz_k_points(self):
1175 raise NotImplementedError
1177 def read_vib_freq(self, lines=None) -> Tuple[List[float], List[float]]:
1178 """Read vibrational frequencies.
1180 Returns:
1181 List of real and list of imaginary frequencies
1182 (imaginary number as real number).
1183 """
1184 freq = []
1185 i_freq = []
1187 if not lines:
1188 lines = self.load_file('OUTCAR')
1190 for line in lines:
1191 data = line.split()
1192 if 'THz' in data:
1193 if 'f/i=' not in data:
1194 freq.append(float(data[-2]))
1195 else:
1196 i_freq.append(float(data[-2]))
1197 return freq, i_freq
1199 def _read_massweighted_hessian_xml(self) -> np.ndarray:
1200 """Read the Mass Weighted Hessian from vasprun.xml.
1202 Returns:
1203 The Mass Weighted Hessian as np.ndarray from the xml file.
1205 Raises a ReadError if the reader is not able to read the Hessian.
1207 Converts to ASE units for VASP version 6.
1208 """
1210 file = self._indir('vasprun.xml')
1211 try:
1212 tree = ElementTree.iterparse(file)
1213 hessian = None
1214 for event, elem in tree:
1215 if elem.tag == 'dynmat':
1216 for i, entry in enumerate(
1217 elem.findall('varray[@name="hessian"]/v')):
1218 text_split = entry.text.split()
1219 if not text_split:
1220 raise ElementTree.ParseError(
1221 "Could not find varray hessian!")
1222 if i == 0:
1223 n_items = len(text_split)
1224 hessian = np.zeros((n_items, n_items))
1225 assert isinstance(hessian, np.ndarray)
1226 hessian[i, :] = np.array(
1227 [float(val) for val in text_split])
1228 if i != n_items - 1:
1229 raise ElementTree.ParseError(
1230 "Hessian is not quadratic!")
1231 # VASP6+ uses THz**2 as unit, not mEV**2 as before
1232 for entry in elem.findall('i[@name="unit"]'):
1233 if entry.text.strip() == 'THz^2':
1234 conv = ase.units._amu / ase.units._e / \
1235 1e-4 * (2 * np.pi)**2 # THz**2 to eV**2
1236 # VASP6 uses factor 2pi
1237 # 1e-4 = (angstrom to meter times Hz to THz) squared
1238 # = (1e10 times 1e-12)**2
1239 break
1240 else: # Catch changes in VASP
1241 vasp_version_error_msg = (
1242 f'The file "{file}" is from a '
1243 'non-supported VASP version. '
1244 'Not sure what unit the Hessian '
1245 'is in, aborting.')
1246 raise calculator.ReadError(vasp_version_error_msg)
1248 else:
1249 conv = 1.0 # VASP version <6 unit is meV**2
1250 assert isinstance(hessian, np.ndarray)
1251 hessian *= conv
1252 if hessian is None:
1253 raise ElementTree.ParseError("Hessian is None!")
1255 except ElementTree.ParseError as exc:
1256 incomplete_msg = (
1257 f'The file "{file}" is incomplete, '
1258 'and no DFT data was available. '
1259 'This is likely due to an incomplete calculation.')
1260 raise calculator.ReadError(incomplete_msg) from exc
1261 # VASP uses the negative definition of the hessian compared to ASE
1262 return -hessian
1264 def get_vibrations(self) -> VibrationsData:
1265 """Get a VibrationsData Object from a VASP Calculation.
1267 Returns:
1268 VibrationsData object.
1270 Note that the atoms in the VibrationsData object can be resorted.
1272 Uses the (mass weighted) Hessian from vasprun.xml, different masses
1273 in the POTCAR can therefore result in different results.
1275 Note the limitations concerning k-points and symmetry mentioned in
1276 the VASP-Wiki.
1277 """
1279 mass_weighted_hessian = self._read_massweighted_hessian_xml()
1280 # get indices of freely moving atoms, i.e. respect constraints.
1281 indices = VibrationsData.indices_from_constraints(self.atoms)
1282 # save the corresponding sorted atom numbers
1283 sort_indices = np.array(self.sort)[indices]
1284 # mass weights = 1/sqrt(mass)
1285 mass_weights = np.repeat(self.atoms.get_masses()[sort_indices]**-0.5, 3)
1286 # get the unweighted hessian = H_w / m_w / m_w^T
1287 # ugly and twice the work, but needed since vasprun.xml does
1288 # not have the unweighted ase.vibrations.vibration will do the
1289 # opposite in Vibrations.read
1290 hessian = mass_weighted_hessian / \
1291 mass_weights / mass_weights[:, np.newaxis]
1293 return VibrationsData.from_2d(self.atoms[self.sort], hessian, indices)
1295 def get_nonselfconsistent_energies(self, bee_type):
1296 """ Method that reads and returns BEE energy contributions
1297 written in OUTCAR file.
1298 """
1299 assert bee_type == 'beefvdw'
1300 cmd = 'grep -32 "BEEF xc energy contributions" OUTCAR | tail -32'
1301 p = os.popen(cmd, 'r')
1302 s = p.readlines()
1303 p.close()
1304 xc = np.array([])
1305 for line in s:
1306 l_ = float(line.split(":")[-1])
1307 xc = np.append(xc, l_)
1308 assert len(xc) == 32
1309 return xc
1312#####################################
1313# Below defines helper functions
1314# for the VASP calculator
1315#####################################
1318def check_atoms(atoms: ase.Atoms) -> None:
1319 """Perform checks on the atoms object, to verify that
1320 it can be run by VASP.
1321 A CalculatorSetupError error is raised if the atoms are not supported.
1322 """
1324 # Loop through all check functions
1325 for check in (check_atoms_type, check_cell, check_pbc):
1326 check(atoms)
1329def check_cell(atoms: ase.Atoms) -> None:
1330 """Check if there is a zero unit cell.
1331 Raises CalculatorSetupError if the cell is wrong.
1332 """
1333 if atoms.cell.rank < 3:
1334 raise calculator.CalculatorSetupError(
1335 "The lattice vectors are zero! "
1336 "This is the default value - please specify a "
1337 "unit cell.")
1340def check_pbc(atoms: ase.Atoms) -> None:
1341 """Check if any boundaries are not PBC, as VASP
1342 cannot handle non-PBC.
1343 Raises CalculatorSetupError.
1344 """
1345 if not atoms.pbc.all():
1346 raise calculator.CalculatorSetupError(
1347 "Vasp cannot handle non-periodic boundaries. "
1348 "Please enable all PBC, e.g. atoms.pbc=True")
1351def check_atoms_type(atoms: ase.Atoms) -> None:
1352 """Check that the passed atoms object is in fact an Atoms object.
1353 Raises CalculatorSetupError.
1354 """
1355 if not isinstance(atoms, ase.Atoms):
1356 raise calculator.CalculatorSetupError(
1357 'Expected an Atoms object, '
1358 'instead got object of type {}'.format(type(atoms)))