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