1# Copyright (C) 2008 CSC - Scientific Computing Ltd.
2"""This module defines an ASE interface to VASP.
3
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.
8
9The user should also set the environmental flag $VASP_SCRIPT pointing
10to a python script looking something like::
11
12   import os
13   exitcode = os.system('vasp')
14
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'
17
18http://cms.mpi.univie.ac.at/vasp/
19"""
20
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
31
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
39
40
41class Vasp(GenerateVaspInput, Calculator):  # type: ignore
42    """ASE interface for the Vienna Ab initio Simulation Package (VASP),
43    with the Calculator interface.
44
45        Parameters:
46
47            atoms:  object
48                Attach an atoms object to the calculator.
49
50            label: str
51                Prefix for the output file, and sets the working directory.
52                Default is 'vasp'.
53
54            directory: str
55                Set the working directory. Is prepended to ``label``.
56
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.
61
62            txt: bool, None, str or writable object
63                - If txt is None, output stream will be supressed
64
65                - If txt is '-' the output will be sent through stdout
66
67                - If txt is a string a file will be opened,\
68                    and the output will be sent to that file.
69
70                - Finally, txt can also be a an output stream,\
71                    which has a 'write' attribute.
72
73                Default is 'vasp.out'
74
75                - Examples:
76
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
81
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
88
89    # Environment commands
90    env_commands = ('ASE_VASP_COMMAND', 'VASP_COMMAND', 'VASP_SCRIPT')
91
92    implemented_properties = [
93        'energy', 'free_energy', 'forces', 'dipole', 'fermi', 'stress',
94        'magmom', 'magmoms'
95    ]
96
97    # Can be used later to set some ASE defaults
98    default_parameters: Dict[str, Any] = {}
99
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):
109
110        self._atoms = None
111        self.results = {}
112
113        # Initialize parameter dictionaries
114        GenerateVaspInput.__init__(self)
115        self._store_param_state()  # Initialize an empty parameter state
116
117        # Store calculator from vasprun.xml here - None => uninitialized
118        self._xml_calc = None
119
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
133
134        if isinstance(restart, bool):
135            if restart is True:
136                restart = self.label
137            else:
138                restart = None
139
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)
148
149        self.command = command
150
151        self._txt = None
152        self.txt = txt  # Set the output txt stream
153        self.version = None
154
155        # XXX: This seems to break restarting, unless we return first.
156        # Do we really still need to enfore this?
157
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})
169
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
193
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.
198
199        Allows for setting ``label``, ``directory`` and ``txt``
200        without resetting the results in the calculator.
201        """
202        changed_parameters = {}
203
204        if 'label' in kwargs:
205            self.label = kwargs.pop('label')
206
207        if 'directory' in kwargs:
208            # str() call to deal with pathlib objects
209            self.directory = str(kwargs.pop('directory'))
210
211        if 'txt' in kwargs:
212            self.txt = kwargs.pop('txt')
213
214        if 'atoms' in kwargs:
215            atoms = kwargs.pop('atoms')
216            self.atoms = atoms  # Resets results
217
218        if 'command' in kwargs:
219            self.command = kwargs.pop('command')
220
221        changed_parameters.update(Calculator.set(self, **kwargs))
222
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()
230
231    def reset(self):
232        self.atoms = None
233        self.clear_results()
234
235    def clear_results(self):
236        self.results.clear()
237        self._xml_calc = None
238
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.
246
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'
252
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()
259
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        """
265
266        txt = self.txt
267        open_and_close = False  # Do we open the file?
268
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))
288
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()
296
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.
302
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)
310
311        self.clear_results()
312
313        if atoms is not None:
314            self.atoms = atoms.copy()
315
316        command = self.make_command(self.command)
317        self.write_input(self.atoms, properties, system_changes)
318
319        with self._txt_outstream() as out:
320            errorcode = self._run(command=command,
321                                  out=out,
322                                  directory=self.directory)
323
324        if errorcode:
325            raise calculator.CalculationFailed(
326                '{} in {} returned an error: {:d}'.format(
327                    self.name, self.directory, errorcode))
328
329        # Read results from calculation
330        self.update_atoms(atoms)
331        self.read_results()
332
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
344
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
353
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
359
360        # First we check for default changes
361        system_changes = Calculator.check_state(self, atoms, tol=tol)
362
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)
369
370        return system_changes
371
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())
385
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:
391
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()
401
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        }
409
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        }
417
418        if self.atoms:
419            # Encode atoms as dict
420            from ase.db.row import atoms2dict
421            dct['atoms'] = atoms2dict(self.atoms)
422
423        return dct
424
425    def fromdict(self, dct):
426        """Restore calculator from a :func:`~ase.calculators.vasp.Vasp.asdict`
427        dictionary.
428
429        Parameters:
430
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'])
445
446    def write_json(self, filename):
447        """Dump calculator state to JSON file.
448
449        Parameters:
450
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)
458
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)
463
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)
470
471        self.initialize(atoms)
472
473        GenerateVaspInput.write_input(self, atoms, directory=self.directory)
474
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)
482
483        # If we restart, self.parameters isn't initialized
484        if self.parameters is None:
485            self.parameters = self.get_default_parameters()
486
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))
493
494        # Build sorting and resorting lists
495        self.read_sort()
496
497        # Read atoms
498        self.atoms = self.read_atoms(filename=self._indir('CONTCAR'))
499
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'))
504
505        # Read the results from the calculation
506        self.read_results()
507
508    def _indir(self, filename):
509        """Prepend current directory to filename"""
510        return Path(self.directory) / filename
511
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)
529
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]
534
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
545
546        self.atoms = atoms  # Creates a copy
547
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')
552
553        # Read the data we can from vasprun.xml
554        calc_xml = self._read_xml()
555        xml_results = calc_xml.results
556
557        # Fix sorting
558        xml_results['forces'] = xml_results['forces'][self.resort]
559
560        self.results.update(xml_results)
561
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)
567
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.
571
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))
579
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))
584
585        self._set_old_keywords()
586
587        # Store the parameters used for this calculation
588        self._store_param_state()
589
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()
601
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']
607
608    @kpts.setter
609    def kpts(self, kpts):
610        """Set kpts in input_params dict"""
611        self.input_params['kpts'] = kpts
612
613    @property
614    def encut(self):
615        """Direct access to the encut parameter"""
616        return self.float_params['encut']
617
618    @encut.setter
619    def encut(self, encut):
620        """Direct access for setting the encut parameter"""
621        self.set(encut=encut)
622
623    @property
624    def xc(self):
625        """Direct access to the xc parameter"""
626        return self.get_xc_functional()
627
628    @xc.setter
629    def xc(self, xc):
630        """Direct access for setting the xc parameter"""
631        self.set(xc=xc)
632
633    @property
634    def atoms(self):
635        return self._atoms
636
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()
646
647    def load_file(self, filename):
648        """Reads a file in the directory, and returns the lines
649
650        Example:
651        >>> outcar = load_file('OUTCAR')
652        """
653        filename = self._indir(filename)
654        with open(filename, 'r') as fd:
655            return fd.readlines()
656
657    @contextmanager
658    def load_file_iter(self, filename):
659        """Return a file iterator"""
660
661        filename = self._indir(filename)
662        with open(filename, 'r') as fd:
663            yield fd
664
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()
672
673        self.version = self.get_version()
674
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)
679
680        self.dipole = self.read_dipole(lines=lines)
681
682        self.stress = self.read_stress(lines=lines)
683        self.nbands = self.read_nbands(lines=lines)
684
685        self.read_ldau()
686        self.magnetic_moment, self.magnetic_moments = self.read_mag(
687            lines=lines)
688
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
704
705        if _xml_atoms is None or _xml_atoms.calc is None:
706            raise calculator.ReadError(incomplete_msg)
707
708        self._xml_calc = _xml_atoms.calc
709        return self._xml_calc
710
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
717
718    @_xml_calc.setter
719    def _xml_calc(self, value):
720        self.__xml_calc = value
721
722    def get_ibz_k_points(self):
723        calc = self._xml_calc
724        return calc.get_ibz_k_points()
725
726    def get_kpt(self, kpt=0, spin=0):
727        calc = self._xml_calc
728        return calc.get_kpt(kpt=kpt, spin=spin)
729
730    def get_eigenvalues(self, kpt=0, spin=0):
731        calc = self._xml_calc
732        return calc.get_eigenvalues(kpt=kpt, spin=spin)
733
734    def get_fermi_level(self):
735        calc = self._xml_calc
736        return calc.get_fermi_level()
737
738    def get_homo_lumo(self):
739        calc = self._xml_calc
740        return calc.get_homo_lumo()
741
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)
745
746    def get_occupation_numbers(self, kpt=0, spin=0):
747        calc = self._xml_calc
748        return calc.get_occupation_numbers(kpt, spin)
749
750    def get_spin_polarized(self):
751        calc = self._xml_calc
752        return calc.get_spin_polarized()
753
754    def get_number_of_spins(self):
755        calc = self._xml_calc
756        return calc.get_number_of_spins()
757
758    def get_number_of_bands(self):
759        return self.results.get('nbands', None)
760
761    def get_number_of_electrons(self, lines=None):
762        if not lines:
763            lines = self.load_file('OUTCAR')
764
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
771
772    def get_k_point_weights(self):
773        filename = self._indir('IBZKPT')
774        return self.read_k_point_weights(filename)
775
776    def get_dos(self, spin=None, **kwargs):
777        """
778        The total DOS.
779
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
788
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
794
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
807
808    def get_number_of_iterations(self):
809        return self.read_number_of_iterations()
810
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
819
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
827
828    def read_stress(self, lines=None):
829        """Read stress from OUTCAR.
830
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')
837
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
844
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')
849
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
878
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
884
885    def get_xc_functional(self):
886        """Returns the XC functional or the pseudopotential type
887
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.')
898
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')
905
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]
924
925    def read_forces(self, all=False, lines=None):
926        """Method that reads forces from OUTCAR file.
927
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."""
931
932        if not lines:
933            lines = self.load_file('OUTCAR')
934
935        if all:
936            all_forces = []
937
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]]))
945
946                if all:
947                    all_forces.append(np.array(forces)[self.resort])
948
949        if all:
950            return np.array(all_forces)
951        return np.array(forces)[self.resort]
952
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')
957
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
963
964    def read_dipole(self, lines=None):
965        """Read dipole from OUTCAR"""
966        if not lines:
967            lines = self.load_file('OUTCAR')
968
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
974
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
994
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')
1000
1001        magnetic_moments = np.zeros(len(self.atoms))
1002        magstr = 'magnetization (x)'
1003
1004        # Search for the last occurrence
1005        nidx = -1
1006        for n, line in enumerate(lines):
1007            if magstr in line:
1008                nidx = n
1009
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]
1015
1016    def _read_magnetic_moment(self, lines=None):
1017        """Read magnetic moment from OUTCAR"""
1018        if not lines:
1019            lines = self.load_file('OUTCAR')
1020
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
1025
1026    def read_nbands(self, lines=None):
1027        """Read number of bands from OUTCAR"""
1028        if not lines:
1029            lines = self.load_file('OUTCAR')
1030
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
1036
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')
1041
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
1083
1084    def read_k_point_weights(self, filename):
1085        """Read k-point weighting. Normally named IBZKPT."""
1086
1087        lines = self.load_file(filename)
1088
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)
1098
1099        return kpt_weights
1100
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
1109
1110    def read_spinpol(self, lines=None):
1111        """Method which reads if a calculation from spinpolarized using OUTCAR.
1112
1113        Depreciated: Use get_spin_polarized() instead.
1114        """
1115        if not lines:
1116            lines = self.load_file('OUTCAR')
1117
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
1125
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
1131
1132    @property
1133    def txt(self):
1134        return self._txt
1135
1136    @txt.setter
1137    def txt(self, txt):
1138        if isinstance(txt, PurePath):
1139            txt = str(txt)
1140        self._txt = txt
1141
1142    def get_number_of_grid_points(self):
1143        raise NotImplementedError
1144
1145    def get_pseudo_density(self):
1146        raise NotImplementedError
1147
1148    def get_pseudo_wavefunction(self, n=0, k=0, s=0, pad=True):
1149        raise NotImplementedError
1150
1151    def get_bz_k_points(self):
1152        raise NotImplementedError
1153
1154    def read_vib_freq(self, lines=None):
1155        """Read vibrational frequencies.
1156
1157        Returns list of real and list of imaginary frequencies."""
1158        freq = []
1159        i_freq = []
1160
1161        if not lines:
1162            lines = self.load_file('OUTCAR')
1163
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
1172
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
1188
1189
1190#####################################
1191# Below defines helper functions
1192# for the VASP calculator
1193#####################################
1194
1195
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    """
1201
1202    # Loop through all check functions
1203    for check in (check_atoms_type, check_cell, check_pbc):
1204        check(atoms)
1205
1206
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.")
1216
1217
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")
1227
1228
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))))
1237