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