Coverage for /builds/debichem-team/python-ase/ase/calculators/vasp/vasp.py : 29.57%
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
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
41class Vasp(GenerateVaspInput, Calculator): # type: ignore
42 """ASE interface for the Vienna Ab initio Simulation Package (VASP),
43 with the Calculator interface.
45 Parameters:
47 atoms: object
48 Attach an atoms object to the calculator.
50 label: str
51 Prefix for the output file, and sets the working directory.
52 Default is 'vasp'.
54 directory: str
55 Set the working directory. Is prepended to ``label``.
57 restart: str or bool
58 Sets a label for the directory to load files from.
59 if :code:`restart=True`, the working directory from
60 ``directory`` is used.
62 txt: bool, None, str or writable object
63 - If txt is None, output stream will be supressed
65 - If txt is '-' the output will be sent through stdout
67 - If txt is a string a file will be opened,\
68 and the output will be sent to that file.
70 - Finally, txt can also be a an output stream,\
71 which has a 'write' attribute.
73 Default is 'vasp.out'
75 - Examples:
77 >>> Vasp(label='mylabel', txt='vasp.out') # Redirect stdout
78 >>> Vasp(txt='myfile.txt') # Redirect stdout
79 >>> Vasp(txt='-') # Print vasp output to stdout
80 >>> Vasp(txt=None) # Suppress txt output
82 command: str
83 Custom instructions on how to execute VASP. Has priority over
84 environment variables.
85 """
86 name = 'vasp'
87 ase_objtype = 'vasp_calculator' # For JSON storage
89 # Environment commands
90 env_commands = ('ASE_VASP_COMMAND', 'VASP_COMMAND', 'VASP_SCRIPT')
92 implemented_properties = [
93 'energy', 'free_energy', 'forces', 'dipole', 'fermi', 'stress',
94 'magmom', 'magmoms'
95 ]
97 # Can be used later to set some ASE defaults
98 default_parameters: Dict[str, Any] = {}
100 def __init__(self,
101 atoms=None,
102 restart=None,
103 directory='.',
104 label='vasp',
105 ignore_bad_restart_file=Calculator._deprecated,
106 command=None,
107 txt='vasp.out',
108 **kwargs):
110 self._atoms = None
111 self.results = {}
113 # Initialize parameter dictionaries
114 GenerateVaspInput.__init__(self)
115 self._store_param_state() # Initialize an empty parameter state
117 # Store calculator from vasprun.xml here - None => uninitialized
118 self._xml_calc = None
120 # Set directory and label
121 self.directory = directory
122 if '/' in label:
123 warn(('Specifying directory in "label" is deprecated, '
124 'use "directory" instead.'), np.VisibleDeprecationWarning)
125 if self.directory != '.':
126 raise ValueError('Directory redundantly specified though '
127 'directory="{}" and label="{}". '
128 'Please omit "/" in label.'.format(
129 self.directory, label))
130 self.label = label
131 else:
132 self.prefix = label # The label should only contain the prefix
134 if isinstance(restart, bool):
135 if restart is True:
136 restart = self.label
137 else:
138 restart = None
140 Calculator.__init__(
141 self,
142 restart=restart,
143 ignore_bad_restart_file=ignore_bad_restart_file,
144 # We already, manually, created the label
145 label=self.label,
146 atoms=atoms,
147 **kwargs)
149 self.command = command
151 self._txt = None
152 self.txt = txt # Set the output txt stream
153 self.version = None
155 # XXX: This seems to break restarting, unless we return first.
156 # Do we really still need to enfore this?
158 # # If no XC combination, GGA functional or POTCAR type is specified,
159 # # default to PW91. This is mostly chosen for backwards compatibility.
160 # if kwargs.get('xc', None):
161 # pass
162 # elif not (kwargs.get('gga', None) or kwargs.get('pp', None)):
163 # self.input_params.update({'xc': 'PW91'})
164 # # A null value of xc is permitted; custom recipes can be
165 # # used by explicitly setting the pseudopotential set and
166 # # INCAR keys
167 # else:
168 # self.input_params.update({'xc': None})
170 def make_command(self, command=None):
171 """Return command if one is passed, otherwise try to find
172 ASE_VASP_COMMAND, VASP_COMMAND or VASP_SCRIPT.
173 If none are set, a CalculatorSetupError is raised"""
174 if command:
175 cmd = command
176 else:
177 # Search for the environment commands
178 for env in self.env_commands:
179 if env in os.environ:
180 cmd = os.environ[env].replace('PREFIX', self.prefix)
181 if env == 'VASP_SCRIPT':
182 # Make the system python exe run $VASP_SCRIPT
183 exe = sys.executable
184 cmd = ' '.join([exe, cmd])
185 break
186 else:
187 msg = ('Please set either command in calculator'
188 ' or one of the following environment '
189 'variables (prioritized as follows): {}').format(
190 ', '.join(self.env_commands))
191 raise calculator.CalculatorSetupError(msg)
192 return cmd
194 def set(self, **kwargs):
195 """Override the set function, to test for changes in the
196 Vasp Calculator, then call the create_input.set()
197 on remaining inputs for VASP specific keys.
199 Allows for setting ``label``, ``directory`` and ``txt``
200 without resetting the results in the calculator.
201 """
202 changed_parameters = {}
204 if 'label' in kwargs:
205 self.label = kwargs.pop('label')
207 if 'directory' in kwargs:
208 # str() call to deal with pathlib objects
209 self.directory = str(kwargs.pop('directory'))
211 if 'txt' in kwargs:
212 self.txt = kwargs.pop('txt')
214 if 'atoms' in kwargs:
215 atoms = kwargs.pop('atoms')
216 self.atoms = atoms # Resets results
218 if 'command' in kwargs:
219 self.command = kwargs.pop('command')
221 changed_parameters.update(Calculator.set(self, **kwargs))
223 # We might at some point add more to changed parameters, or use it
224 if changed_parameters:
225 self.clear_results() # We don't want to clear atoms
226 if kwargs:
227 # If we make any changes to Vasp input, we always reset
228 GenerateVaspInput.set(self, **kwargs)
229 self.results.clear()
231 def reset(self):
232 self.atoms = None
233 self.clear_results()
235 def clear_results(self):
236 self.results.clear()
237 self._xml_calc = None
239 @contextmanager
240 def _txt_outstream(self):
241 """Custom function for opening a text output stream. Uses self.txt
242 to determine the output stream, and accepts a string or an open
243 writable object.
244 If a string is used, a new stream is opened, and automatically closes
245 the new stream again when exiting.
247 Examples:
248 # Pass a string
249 calc.txt = 'vasp.out'
250 with calc.txt_outstream() as out:
251 calc.run(out=out) # Redirects the stdout to 'vasp.out'
253 # Use an existing stream
254 mystream = open('vasp.out', 'w')
255 calc.txt = mystream
256 with calc.txt_outstream() as out:
257 calc.run(out=out)
258 mystream.close()
260 # Print to stdout
261 calc.txt = '-'
262 with calc.txt_outstream() as out:
263 calc.run(out=out) # output is written to stdout
264 """
266 txt = self.txt
267 open_and_close = False # Do we open the file?
269 if txt is None:
270 # Suppress stdout
271 out = subprocess.DEVNULL
272 else:
273 if isinstance(txt, str):
274 if txt == '-':
275 # subprocess.call redirects this to stdout
276 out = None
277 else:
278 # Open the file in the work directory
279 txt = self._indir(txt)
280 # We wait with opening the file, until we are inside the
281 # try/finally
282 open_and_close = True
283 elif hasattr(txt, 'write'):
284 out = txt
285 else:
286 raise RuntimeError('txt should either be a string'
287 'or an I/O stream, got {}'.format(txt))
289 try:
290 if open_and_close:
291 out = open(txt, 'w')
292 yield out
293 finally:
294 if open_and_close:
295 out.close()
297 def calculate(self,
298 atoms=None,
299 properties=('energy', ),
300 system_changes=tuple(calculator.all_changes)):
301 """Do a VASP calculation in the specified directory.
303 This will generate the necessary VASP input files, and then
304 execute VASP. After execution, the energy, forces. etc. are read
305 from the VASP output files.
306 """
307 # Check for zero-length lattice vectors and PBC
308 # and that we actually have an Atoms object.
309 check_atoms(atoms)
311 self.clear_results()
313 if atoms is not None:
314 self.atoms = atoms.copy()
316 command = self.make_command(self.command)
317 self.write_input(self.atoms, properties, system_changes)
319 with self._txt_outstream() as out:
320 errorcode = self._run(command=command,
321 out=out,
322 directory=self.directory)
324 if errorcode:
325 raise calculator.CalculationFailed(
326 '{} in {} returned an error: {:d}'.format(
327 self.name, self.directory, errorcode))
329 # Read results from calculation
330 self.update_atoms(atoms)
331 self.read_results()
333 def _run(self, command=None, out=None, directory=None):
334 """Method to explicitly execute VASP"""
335 if command is None:
336 command = self.command
337 if directory is None:
338 directory = self.directory
339 errorcode = subprocess.call(command,
340 shell=True,
341 stdout=out,
342 cwd=directory)
343 return errorcode
345 def check_state(self, atoms, tol=1e-15):
346 """Check for system changes since last calculation."""
347 def compare_dict(d1, d2):
348 """Helper function to compare dictionaries"""
349 # Use symmetric difference to find keys which aren't shared
350 # for python 2.7 compatibility
351 if set(d1.keys()) ^ set(d2.keys()):
352 return False
354 # Check for differences in values
355 for key, value in d1.items():
356 if np.any(value != d2[key]):
357 return False
358 return True
360 # First we check for default changes
361 system_changes = Calculator.check_state(self, atoms, tol=tol)
363 # We now check if we have made any changes to the input parameters
364 # XXX: Should we add these parameters to all_changes?
365 for param_string, old_dict in self.param_state.items():
366 param_dict = getattr(self, param_string) # Get current param dict
367 if not compare_dict(param_dict, old_dict):
368 system_changes.append(param_string)
370 return system_changes
372 def _store_param_state(self):
373 """Store current parameter state"""
374 self.param_state = dict(
375 float_params=self.float_params.copy(),
376 exp_params=self.exp_params.copy(),
377 string_params=self.string_params.copy(),
378 int_params=self.int_params.copy(),
379 input_params=self.input_params.copy(),
380 bool_params=self.bool_params.copy(),
381 list_int_params=self.list_int_params.copy(),
382 list_bool_params=self.list_bool_params.copy(),
383 list_float_params=self.list_float_params.copy(),
384 dict_params=self.dict_params.copy())
386 def asdict(self):
387 """Return a dictionary representation of the calculator state.
388 Does NOT contain information on the ``command``, ``txt`` or
389 ``directory`` keywords.
390 Contains the following keys:
392 - ``ase_version``
393 - ``vasp_version``
394 - ``inputs``
395 - ``results``
396 - ``atoms`` (Only if the calculator has an ``Atoms`` object)
397 """
398 # Get versions
399 asevers = ase.__version__
400 vaspvers = self.get_version()
402 self._store_param_state() # Update param state
403 # Store input parameters which have been set
404 inputs = {
405 key: value
406 for param_dct in self.param_state.values()
407 for key, value in param_dct.items() if value is not None
408 }
410 dct = {
411 'ase_version': asevers,
412 'vasp_version': vaspvers,
413 # '__ase_objtype__': self.ase_objtype,
414 'inputs': inputs,
415 'results': self.results.copy()
416 }
418 if self.atoms:
419 # Encode atoms as dict
420 from ase.db.row import atoms2dict
421 dct['atoms'] = atoms2dict(self.atoms)
423 return dct
425 def fromdict(self, dct):
426 """Restore calculator from a :func:`~ase.calculators.vasp.Vasp.asdict`
427 dictionary.
429 Parameters:
431 dct: Dictionary
432 The dictionary which is used to restore the calculator state.
433 """
434 if 'vasp_version' in dct:
435 self.version = dct['vasp_version']
436 if 'inputs' in dct:
437 self.set(**dct['inputs'])
438 self._store_param_state()
439 if 'atoms' in dct:
440 from ase.db.row import AtomsRow
441 atoms = AtomsRow(dct['atoms']).toatoms()
442 self.atoms = atoms
443 if 'results' in dct:
444 self.results.update(dct['results'])
446 def write_json(self, filename):
447 """Dump calculator state to JSON file.
449 Parameters:
451 filename: string
452 The filename which the JSON file will be stored to.
453 Prepends the ``directory`` path to the filename.
454 """
455 filename = self._indir(filename)
456 dct = self.asdict()
457 jsonio.write_json(filename, dct)
459 def read_json(self, filename):
460 """Load Calculator state from an exported JSON Vasp file."""
461 dct = jsonio.read_json(filename)
462 self.fromdict(dct)
464 def write_input(self, atoms, properties=None, system_changes=None):
465 """Write VASP inputfiles, INCAR, KPOINTS and POTCAR"""
466 # Create the folders where we write the files, if we aren't in the
467 # current working directory.
468 if self.directory != os.curdir and not os.path.isdir(self.directory):
469 os.makedirs(self.directory)
471 self.initialize(atoms)
473 GenerateVaspInput.write_input(self, atoms, directory=self.directory)
475 def read(self, label=None):
476 """Read results from VASP output files.
477 Files which are read: OUTCAR, CONTCAR and vasprun.xml
478 Raises ReadError if they are not found"""
479 if label is None:
480 label = self.label
481 Calculator.read(self, label)
483 # If we restart, self.parameters isn't initialized
484 if self.parameters is None:
485 self.parameters = self.get_default_parameters()
487 # Check for existence of the necessary output files
488 for f in ['OUTCAR', 'CONTCAR', 'vasprun.xml']:
489 file = self._indir(f)
490 if not file.is_file():
491 raise calculator.ReadError(
492 'VASP outputfile {} was not found'.format(file))
494 # Build sorting and resorting lists
495 self.read_sort()
497 # Read atoms
498 self.atoms = self.read_atoms(filename=self._indir('CONTCAR'))
500 # Read parameters
501 self.read_incar(filename=self._indir('INCAR'))
502 self.read_kpoints(filename=self._indir('KPOINTS'))
503 self.read_potcar(filename=self._indir('POTCAR'))
505 # Read the results from the calculation
506 self.read_results()
508 def _indir(self, filename):
509 """Prepend current directory to filename"""
510 return Path(self.directory) / filename
512 def read_sort(self):
513 """Create the sorting and resorting list from ase-sort.dat.
514 If the ase-sort.dat file does not exist, the sorting is redone.
515 """
516 sortfile = self._indir('ase-sort.dat')
517 if os.path.isfile(sortfile):
518 self.sort = []
519 self.resort = []
520 with open(sortfile, 'r') as fd:
521 for line in fd:
522 sort, resort = line.split()
523 self.sort.append(int(sort))
524 self.resort.append(int(resort))
525 else:
526 # Redo the sorting
527 atoms = read(self._indir('CONTCAR'))
528 self.initialize(atoms)
530 def read_atoms(self, filename):
531 """Read the atoms from file located in the VASP
532 working directory. Normally called CONTCAR."""
533 return read(filename)[self.resort]
535 def update_atoms(self, atoms):
536 """Update the atoms object with new positions and cell"""
537 if (self.int_params['ibrion'] is not None
538 and self.int_params['nsw'] is not None):
539 if self.int_params['ibrion'] > -1 and self.int_params['nsw'] > 0:
540 # Update atomic positions and unit cell with the ones read
541 # from CONTCAR.
542 atoms_sorted = read(self._indir('CONTCAR'))
543 atoms.positions = atoms_sorted[self.resort].positions
544 atoms.cell = atoms_sorted.cell
546 self.atoms = atoms # Creates a copy
548 def read_results(self):
549 """Read the results from VASP output files"""
550 # Temporarily load OUTCAR into memory
551 outcar = self.load_file('OUTCAR')
553 # Read the data we can from vasprun.xml
554 calc_xml = self._read_xml()
555 xml_results = calc_xml.results
557 # Fix sorting
558 xml_results['forces'] = xml_results['forces'][self.resort]
560 self.results.update(xml_results)
562 # Parse the outcar, as some properties are not loaded in vasprun.xml
563 # We want to limit this as much as possible, as reading large OUTCAR's
564 # is relatively slow
565 # Removed for now
566 # self.read_outcar(lines=outcar)
568 # Update results dict with results from OUTCAR
569 # which aren't written to the atoms object we read from
570 # the vasprun.xml file.
572 self.converged = self.read_convergence(lines=outcar)
573 self.version = self.read_version()
574 magmom, magmoms = self.read_mag(lines=outcar)
575 dipole = self.read_dipole(lines=outcar)
576 nbands = self.read_nbands(lines=outcar)
577 self.results.update(
578 dict(magmom=magmom, magmoms=magmoms, dipole=dipole, nbands=nbands))
580 # Stress is not always present.
581 # Prevent calculation from going into a loop
582 if 'stress' not in self.results:
583 self.results.update(dict(stress=None))
585 self._set_old_keywords()
587 # Store the parameters used for this calculation
588 self._store_param_state()
590 def _set_old_keywords(self):
591 """Store keywords for backwards compatibility wd VASP calculator"""
592 self.spinpol = self.get_spin_polarized()
593 self.energy_free = self.get_potential_energy(force_consistent=True)
594 self.energy_zero = self.get_potential_energy(force_consistent=False)
595 self.forces = self.get_forces()
596 self.fermi = self.get_fermi_level()
597 self.dipole = self.get_dipole_moment()
598 # Prevent calculation from going into a loop
599 self.stress = self.get_property('stress', allow_calculation=False)
600 self.nbands = self.get_number_of_bands()
602 # Below defines some functions for faster access to certain common keywords
603 @property
604 def kpts(self):
605 """Access the kpts from input_params dict"""
606 return self.input_params['kpts']
608 @kpts.setter
609 def kpts(self, kpts):
610 """Set kpts in input_params dict"""
611 self.input_params['kpts'] = kpts
613 @property
614 def encut(self):
615 """Direct access to the encut parameter"""
616 return self.float_params['encut']
618 @encut.setter
619 def encut(self, encut):
620 """Direct access for setting the encut parameter"""
621 self.set(encut=encut)
623 @property
624 def xc(self):
625 """Direct access to the xc parameter"""
626 return self.get_xc_functional()
628 @xc.setter
629 def xc(self, xc):
630 """Direct access for setting the xc parameter"""
631 self.set(xc=xc)
633 @property
634 def atoms(self):
635 return self._atoms
637 @atoms.setter
638 def atoms(self, atoms):
639 if atoms is None:
640 self._atoms = None
641 self.clear_results()
642 else:
643 if self.check_state(atoms):
644 self.clear_results()
645 self._atoms = atoms.copy()
647 def load_file(self, filename):
648 """Reads a file in the directory, and returns the lines
650 Example:
651 >>> outcar = load_file('OUTCAR')
652 """
653 filename = self._indir(filename)
654 with open(filename, 'r') as fd:
655 return fd.readlines()
657 @contextmanager
658 def load_file_iter(self, filename):
659 """Return a file iterator"""
661 filename = self._indir(filename)
662 with open(filename, 'r') as fd:
663 yield fd
665 def read_outcar(self, lines=None):
666 """Read results from the OUTCAR file.
667 Deprecated, see read_results()"""
668 if not lines:
669 lines = self.load_file('OUTCAR')
670 # Spin polarized calculation?
671 self.spinpol = self.get_spin_polarized()
673 self.version = self.get_version()
675 # XXX: Do we want to read all of this again?
676 self.energy_free, self.energy_zero = self.read_energy(lines=lines)
677 self.forces = self.read_forces(lines=lines)
678 self.fermi = self.read_fermi(lines=lines)
680 self.dipole = self.read_dipole(lines=lines)
682 self.stress = self.read_stress(lines=lines)
683 self.nbands = self.read_nbands(lines=lines)
685 self.read_ldau()
686 self.magnetic_moment, self.magnetic_moments = self.read_mag(
687 lines=lines)
689 def _read_xml(self) -> SinglePointDFTCalculator:
690 """Read vasprun.xml, and return the last calculator object.
691 Returns calculator from the xml file.
692 Raises a ReadError if the reader is not able to construct a calculator.
693 """
694 file = self._indir('vasprun.xml')
695 incomplete_msg = (
696 f'The file "{file}" is incomplete, and no DFT data was available. '
697 'This is likely due to an incomplete calculation.')
698 try:
699 _xml_atoms = read(file, index=-1, format='vasp-xml')
700 # Silence mypy, we should only ever get a single atoms object
701 assert isinstance(_xml_atoms, ase.Atoms)
702 except ElementTree.ParseError as exc:
703 raise calculator.ReadError(incomplete_msg) from exc
705 if _xml_atoms is None or _xml_atoms.calc is None:
706 raise calculator.ReadError(incomplete_msg)
708 self._xml_calc = _xml_atoms.calc
709 return self._xml_calc
711 @property
712 def _xml_calc(self) -> SinglePointDFTCalculator:
713 if self.__xml_calc is None:
714 raise RuntimeError(('vasprun.xml data has not yet been loaded. '
715 'Run read_results() first.'))
716 return self.__xml_calc
718 @_xml_calc.setter
719 def _xml_calc(self, value):
720 self.__xml_calc = value
722 def get_ibz_k_points(self):
723 calc = self._xml_calc
724 return calc.get_ibz_k_points()
726 def get_kpt(self, kpt=0, spin=0):
727 calc = self._xml_calc
728 return calc.get_kpt(kpt=kpt, spin=spin)
730 def get_eigenvalues(self, kpt=0, spin=0):
731 calc = self._xml_calc
732 return calc.get_eigenvalues(kpt=kpt, spin=spin)
734 def get_fermi_level(self):
735 calc = self._xml_calc
736 return calc.get_fermi_level()
738 def get_homo_lumo(self):
739 calc = self._xml_calc
740 return calc.get_homo_lumo()
742 def get_homo_lumo_by_spin(self, spin=0):
743 calc = self._xml_calc
744 return calc.get_homo_lumo_by_spin(spin=spin)
746 def get_occupation_numbers(self, kpt=0, spin=0):
747 calc = self._xml_calc
748 return calc.get_occupation_numbers(kpt, spin)
750 def get_spin_polarized(self):
751 calc = self._xml_calc
752 return calc.get_spin_polarized()
754 def get_number_of_spins(self):
755 calc = self._xml_calc
756 return calc.get_number_of_spins()
758 def get_number_of_bands(self):
759 return self.results.get('nbands', None)
761 def get_number_of_electrons(self, lines=None):
762 if not lines:
763 lines = self.load_file('OUTCAR')
765 nelect = None
766 for line in lines:
767 if 'total number of electrons' in line:
768 nelect = float(line.split('=')[1].split()[0].strip())
769 break
770 return nelect
772 def get_k_point_weights(self):
773 filename = self._indir('IBZKPT')
774 return self.read_k_point_weights(filename)
776 def get_dos(self, spin=None, **kwargs):
777 """
778 The total DOS.
780 Uses the ASE DOS module, and returns a tuple with
781 (energies, dos).
782 """
783 from ase.dft.dos import DOS
784 dos = DOS(self, **kwargs)
785 e = dos.get_energies()
786 d = dos.get_dos(spin=spin)
787 return e, d
789 def get_version(self):
790 if self.version is None:
791 # Try if we can read the version number
792 self.version = self.read_version()
793 return self.version
795 def read_version(self):
796 """Get the VASP version number"""
797 # The version number is the first occurrence, so we can just
798 # load the OUTCAR, as we will return soon anyway
799 if not os.path.isfile(self._indir('OUTCAR')):
800 return None
801 with self.load_file_iter('OUTCAR') as lines:
802 for line in lines:
803 if ' vasp.' in line:
804 return line[len(' vasp.'):].split()[0]
805 # We didn't find the version in VASP
806 return None
808 def get_number_of_iterations(self):
809 return self.read_number_of_iterations()
811 def read_number_of_iterations(self):
812 niter = None
813 with self.load_file_iter('OUTCAR') as lines:
814 for line in lines:
815 # find the last iteration number
816 if '- Iteration' in line:
817 niter = list(map(int, re.findall(r'\d+', line)))[1]
818 return niter
820 def read_number_of_ionic_steps(self):
821 niter = None
822 with self.load_file_iter('OUTCAR') as lines:
823 for line in lines:
824 if '- Iteration' in line:
825 niter = list(map(int, re.findall(r'\d+', line)))[0]
826 return niter
828 def read_stress(self, lines=None):
829 """Read stress from OUTCAR.
831 Depreciated: Use get_stress() instead.
832 """
833 # We don't really need this, as we read this from vasprun.xml
834 # keeping it around "just in case" for now
835 if not lines:
836 lines = self.load_file('OUTCAR')
838 stress = None
839 for line in lines:
840 if ' in kB ' in line:
841 stress = -np.array([float(a) for a in line.split()[2:]])
842 stress = stress[[0, 1, 2, 4, 5, 3]] * 1e-1 * ase.units.GPa
843 return stress
845 def read_ldau(self, lines=None):
846 """Read the LDA+U values from OUTCAR"""
847 if not lines:
848 lines = self.load_file('OUTCAR')
850 ldau_luj = None
851 ldauprint = None
852 ldau = None
853 ldautype = None
854 atomtypes = []
855 # read ldau parameters from outcar
856 for line in lines:
857 if line.find('TITEL') != -1: # What atoms are present
858 atomtypes.append(line.split()[3].split('_')[0].split('.')[0])
859 if line.find('LDAUTYPE') != -1: # Is this a DFT+U calculation
860 ldautype = int(line.split('=')[-1])
861 ldau = True
862 ldau_luj = {}
863 if line.find('LDAUL') != -1:
864 L = line.split('=')[-1].split()
865 if line.find('LDAUU') != -1:
866 U = line.split('=')[-1].split()
867 if line.find('LDAUJ') != -1:
868 J = line.split('=')[-1].split()
869 # create dictionary
870 if ldau:
871 for i, symbol in enumerate(atomtypes):
872 ldau_luj[symbol] = {
873 'L': int(L[i]),
874 'U': float(U[i]),
875 'J': float(J[i])
876 }
877 self.dict_params['ldau_luj'] = ldau_luj
879 self.ldau = ldau
880 self.ldauprint = ldauprint
881 self.ldautype = ldautype
882 self.ldau_luj = ldau_luj
883 return ldau, ldauprint, ldautype, ldau_luj
885 def get_xc_functional(self):
886 """Returns the XC functional or the pseudopotential type
888 If a XC recipe is set explicitly with 'xc', this is returned.
889 Otherwise, the XC functional associated with the
890 pseudopotentials (LDA, PW91 or PBE) is returned.
891 The string is always cast to uppercase for consistency
892 in checks."""
893 if self.input_params.get('xc', None):
894 return self.input_params['xc'].upper()
895 if self.input_params.get('pp', None):
896 return self.input_params['pp'].upper()
897 raise ValueError('No xc or pp found.')
899 # Methods for reading information from OUTCAR files:
900 def read_energy(self, all=None, lines=None):
901 """Method to read energy from OUTCAR file.
902 Depreciated: use get_potential_energy() instead"""
903 if not lines:
904 lines = self.load_file('OUTCAR')
906 [energy_free, energy_zero] = [0, 0]
907 if all:
908 energy_free = []
909 energy_zero = []
910 for line in lines:
911 # Free energy
912 if line.lower().startswith(' free energy toten'):
913 if all:
914 energy_free.append(float(line.split()[-2]))
915 else:
916 energy_free = float(line.split()[-2])
917 # Extrapolated zero point energy
918 if line.startswith(' energy without entropy'):
919 if all:
920 energy_zero.append(float(line.split()[-1]))
921 else:
922 energy_zero = float(line.split()[-1])
923 return [energy_free, energy_zero]
925 def read_forces(self, all=False, lines=None):
926 """Method that reads forces from OUTCAR file.
928 If 'all' is switched on, the forces for all ionic steps
929 in the OUTCAR file be returned, in other case only the
930 forces for the last ionic configuration is returned."""
932 if not lines:
933 lines = self.load_file('OUTCAR')
935 if all:
936 all_forces = []
938 for n, line in enumerate(lines):
939 if 'TOTAL-FORCE' in line:
940 forces = []
941 for i in range(len(self.atoms)):
942 forces.append(
943 np.array(
944 [float(f) for f in lines[n + 2 + i].split()[3:6]]))
946 if all:
947 all_forces.append(np.array(forces)[self.resort])
949 if all:
950 return np.array(all_forces)
951 return np.array(forces)[self.resort]
953 def read_fermi(self, lines=None):
954 """Method that reads Fermi energy from OUTCAR file"""
955 if not lines:
956 lines = self.load_file('OUTCAR')
958 E_f = None
959 for line in lines:
960 if 'E-fermi' in line:
961 E_f = float(line.split()[2])
962 return E_f
964 def read_dipole(self, lines=None):
965 """Read dipole from OUTCAR"""
966 if not lines:
967 lines = self.load_file('OUTCAR')
969 dipolemoment = np.zeros([1, 3])
970 for line in lines:
971 if 'dipolmoment' in line:
972 dipolemoment = np.array([float(f) for f in line.split()[1:4]])
973 return dipolemoment
975 def read_mag(self, lines=None):
976 if not lines:
977 lines = self.load_file('OUTCAR')
978 p = self.int_params
979 q = self.list_float_params
980 if self.spinpol:
981 magnetic_moment = self._read_magnetic_moment(lines=lines)
982 if ((p['lorbit'] is not None and p['lorbit'] >= 10)
983 or (p['lorbit'] is None and q['rwigs'])):
984 magnetic_moments = self._read_magnetic_moments(lines=lines)
985 else:
986 warn(('Magnetic moment data not written in OUTCAR (LORBIT<10),'
987 ' setting magnetic_moments to zero.\nSet LORBIT>=10'
988 ' to get information on magnetic moments'))
989 magnetic_moments = np.zeros(len(self.atoms))
990 else:
991 magnetic_moment = 0.0
992 magnetic_moments = np.zeros(len(self.atoms))
993 return magnetic_moment, magnetic_moments
995 def _read_magnetic_moments(self, lines=None):
996 """Read magnetic moments from OUTCAR.
997 Only reads the last occurrence. """
998 if not lines:
999 lines = self.load_file('OUTCAR')
1001 magnetic_moments = np.zeros(len(self.atoms))
1002 magstr = 'magnetization (x)'
1004 # Search for the last occurrence
1005 nidx = -1
1006 for n, line in enumerate(lines):
1007 if magstr in line:
1008 nidx = n
1010 # Read that occurrence
1011 if nidx > -1:
1012 for m in range(len(self.atoms)):
1013 magnetic_moments[m] = float(lines[nidx + m + 4].split()[4])
1014 return magnetic_moments[self.resort]
1016 def _read_magnetic_moment(self, lines=None):
1017 """Read magnetic moment from OUTCAR"""
1018 if not lines:
1019 lines = self.load_file('OUTCAR')
1021 for n, line in enumerate(lines):
1022 if 'number of electron ' in line:
1023 magnetic_moment = float(line.split()[-1])
1024 return magnetic_moment
1026 def read_nbands(self, lines=None):
1027 """Read number of bands from OUTCAR"""
1028 if not lines:
1029 lines = self.load_file('OUTCAR')
1031 for line in lines:
1032 line = self.strip_warnings(line)
1033 if 'NBANDS' in line:
1034 return int(line.split()[-1])
1035 return None
1037 def read_convergence(self, lines=None):
1038 """Method that checks whether a calculation has converged."""
1039 if not lines:
1040 lines = self.load_file('OUTCAR')
1042 converged = None
1043 # First check electronic convergence
1044 for line in lines:
1045 if 0: # vasp always prints that!
1046 if line.rfind('aborting loop') > -1: # scf failed
1047 raise RuntimeError(line.strip())
1048 break
1049 if 'EDIFF ' in line:
1050 ediff = float(line.split()[2])
1051 if 'total energy-change' in line:
1052 # I saw this in an atomic oxygen calculation. it
1053 # breaks this code, so I am checking for it here.
1054 if 'MIXING' in line:
1055 continue
1056 split = line.split(':')
1057 a = float(split[1].split('(')[0])
1058 b = split[1].split('(')[1][0:-2]
1059 # sometimes this line looks like (second number wrong format!):
1060 # energy-change (2. order) :-0.2141803E-08 ( 0.2737684-111)
1061 # we are checking still the first number so
1062 # let's "fix" the format for the second one
1063 if 'e' not in b.lower():
1064 # replace last occurrence of - (assumed exponent) with -e
1065 bsplit = b.split('-')
1066 bsplit[-1] = 'e' + bsplit[-1]
1067 b = '-'.join(bsplit).replace('-e', 'e-')
1068 b = float(b)
1069 if [abs(a), abs(b)] < [ediff, ediff]:
1070 converged = True
1071 else:
1072 converged = False
1073 continue
1074 # Then if ibrion in [1,2,3] check whether ionic relaxation
1075 # condition been fulfilled
1076 if ((self.int_params['ibrion'] in [1, 2, 3]
1077 and self.int_params['nsw'] not in [0])):
1078 if not self.read_relaxed():
1079 converged = False
1080 else:
1081 converged = True
1082 return converged
1084 def read_k_point_weights(self, filename):
1085 """Read k-point weighting. Normally named IBZKPT."""
1087 lines = self.load_file(filename)
1089 if 'Tetrahedra\n' in lines:
1090 N = lines.index('Tetrahedra\n')
1091 else:
1092 N = len(lines)
1093 kpt_weights = []
1094 for n in range(3, N):
1095 kpt_weights.append(float(lines[n].split()[3]))
1096 kpt_weights = np.array(kpt_weights)
1097 kpt_weights /= np.sum(kpt_weights)
1099 return kpt_weights
1101 def read_relaxed(self, lines=None):
1102 """Check if ionic relaxation completed"""
1103 if not lines:
1104 lines = self.load_file('OUTCAR')
1105 for line in lines:
1106 if 'reached required accuracy' in line:
1107 return True
1108 return False
1110 def read_spinpol(self, lines=None):
1111 """Method which reads if a calculation from spinpolarized using OUTCAR.
1113 Depreciated: Use get_spin_polarized() instead.
1114 """
1115 if not lines:
1116 lines = self.load_file('OUTCAR')
1118 for line in lines:
1119 if 'ISPIN' in line:
1120 if int(line.split()[2]) == 2:
1121 self.spinpol = True
1122 else:
1123 self.spinpol = False
1124 return self.spinpol
1126 def strip_warnings(self, line):
1127 """Returns empty string instead of line from warnings in OUTCAR."""
1128 if line[0] == "|":
1129 return ""
1130 return line
1132 @property
1133 def txt(self):
1134 return self._txt
1136 @txt.setter
1137 def txt(self, txt):
1138 if isinstance(txt, PurePath):
1139 txt = str(txt)
1140 self._txt = txt
1142 def get_number_of_grid_points(self):
1143 raise NotImplementedError
1145 def get_pseudo_density(self):
1146 raise NotImplementedError
1148 def get_pseudo_wavefunction(self, n=0, k=0, s=0, pad=True):
1149 raise NotImplementedError
1151 def get_bz_k_points(self):
1152 raise NotImplementedError
1154 def read_vib_freq(self, lines=None):
1155 """Read vibrational frequencies.
1157 Returns list of real and list of imaginary frequencies."""
1158 freq = []
1159 i_freq = []
1161 if not lines:
1162 lines = self.load_file('OUTCAR')
1164 for line in lines:
1165 data = line.split()
1166 if 'THz' in data:
1167 if 'f/i=' not in data:
1168 freq.append(float(data[-2]))
1169 else:
1170 i_freq.append(float(data[-2]))
1171 return freq, i_freq
1173 def get_nonselfconsistent_energies(self, bee_type):
1174 """ Method that reads and returns BEE energy contributions
1175 written in OUTCAR file.
1176 """
1177 assert bee_type == 'beefvdw'
1178 cmd = 'grep -32 "BEEF xc energy contributions" OUTCAR | tail -32'
1179 p = os.popen(cmd, 'r')
1180 s = p.readlines()
1181 p.close()
1182 xc = np.array([])
1183 for line in s:
1184 l_ = float(line.split(":")[-1])
1185 xc = np.append(xc, l_)
1186 assert len(xc) == 32
1187 return xc
1190#####################################
1191# Below defines helper functions
1192# for the VASP calculator
1193#####################################
1196def check_atoms(atoms: ase.Atoms) -> None:
1197 """Perform checks on the atoms object, to verify that
1198 it can be run by VASP.
1199 A CalculatorSetupError error is raised if the atoms are not supported.
1200 """
1202 # Loop through all check functions
1203 for check in (check_atoms_type, check_cell, check_pbc):
1204 check(atoms)
1207def check_cell(atoms: ase.Atoms) -> None:
1208 """Check if there is a zero unit cell.
1209 Raises CalculatorSetupError if the cell is wrong.
1210 """
1211 if atoms.cell.rank < 3:
1212 raise calculator.CalculatorSetupError(
1213 "The lattice vectors are zero! "
1214 "This is the default value - please specify a "
1215 "unit cell.")
1218def check_pbc(atoms: ase.Atoms) -> None:
1219 """Check if any boundaries are not PBC, as VASP
1220 cannot handle non-PBC.
1221 Raises CalculatorSetupError.
1222 """
1223 if not atoms.pbc.all():
1224 raise calculator.CalculatorSetupError(
1225 "Vasp cannot handle non-periodic boundaries. "
1226 "Please enable all PBC, e.g. atoms.pbc=True")
1229def check_atoms_type(atoms: ase.Atoms) -> None:
1230 """Check that the passed atoms object is in fact an Atoms object.
1231 Raises CalculatorSetupError.
1232 """
1233 if not isinstance(atoms, ase.Atoms):
1234 raise calculator.CalculatorSetupError(
1235 ('Expected an Atoms object, '
1236 'instead got object of type {}'.format(type(atoms))))