Merge pull request #1627 from mturiansky/eigenval

Eigenval object

pymatgen/io/vasp/outputs.py

@@ -4499,6 +4499,126 @@ def get_parchg(self, poscar, kpoint, band, spin=None, phase=False,
         return Chgcar(poscar, data)
 
 
+class Eigenval:
+    """
+    Object for reading EIGENVAL file.
+
+    .. attribute:: filename
+
+        string containing input filename
+
+    .. attribute:: occu_tol
+
+        tolerance for determining occupation in band properties
+
+    .. attribute:: ispin
+
+        spin polarization tag (int)
+
+    .. attribute:: nelect
+
+        number of electrons
+
+    .. attribute:: nkpt
+
+        number of kpoints
+
+    .. attribute:: nbands
+
+        number of bands
+
+    .. attribute:: kpoints
+
+        list of kpoints
+
+    .. attribute:: kpoints_weights
+
+        weights of each kpoint in the BZ, should sum to 1.
+
+    .. attribute:: eigenvalues
+
+        Eigenvalues as a dict of {(spin): np.ndarray(shape=(nkpt, nbands, 2))}.
+        This representation is based on actual ordering in VASP and is meant as
+        an intermediate representation to be converted into proper objects. The
+        kpoint index is 0-based (unlike the 1-based indexing in VASP).
+    """
+
+    def __init__(self, filename, occu_tol=1e-8):
+        """
+        Reads input from filename to construct Eigenval object
+
+        Args:
+            filename (str):     filename of EIGENVAL to read in
+            occu_tol (float):   tolerance for determining band gap
+
+        Returns:
+            a pymatgen.io.vasp.outputs.Eigenval object
+        """
+
+        self.filename = filename
+        self.occu_tol = occu_tol
+
+        with zopen(filename, 'r') as f:
+            self.ispin = int(f.readline().split()[-1])
+
+            # useless header information
+            for _ in range(4):
+                f.readline()
+
+            self.nelect, self.nkpt, self.nbands = \
+                list(map(int, f.readline().split()))
+
+            self.kpoints = []
+            self.kpoints_weights = []
+            if self.ispin == 2:
+                self.eigenvalues = \
+                    {Spin.up: np.zeros((self.nkpt, self.nbands, 2)),
+                     Spin.down: np.zeros((self.nkpt, self.nbands, 2))}
+            else:
+                self.eigenvalues = \
+                    {Spin.up: np.zeros((self.nkpt, self.nbands, 2))}
+
+            ikpt = -1
+            for line in f:
+                if re.search(r'(\s+[\-+0-9eE.]+){4}', str(line)):
+                    ikpt += 1
+                    kpt = list(map(float, line.split()))
+                    self.kpoints.append(kpt[:-1])
+                    self.kpoints_weights.append(kpt[-1])
+                    for i in range(self.nbands):
+                        sl = list(map(float, f.readline().split()))
+                        if len(sl) == 3:
+                            self.eigenvalues[Spin.up][ikpt, i, 0] = sl[1]
+                            self.eigenvalues[Spin.up][ikpt, i, 1] = sl[2]
+                        elif len(sl) == 5:
+                            self.eigenvalues[Spin.up][ikpt, i, 0] = sl[1]
+                            self.eigenvalues[Spin.up][ikpt, i, 1] = sl[3]
+                            self.eigenvalues[Spin.down][ikpt, i, 0] = sl[2]
+                            self.eigenvalues[Spin.down][ikpt, i, 1] = sl[4]
+
+    @property
+    def eigenvalue_band_properties(self):
+        """
+        Band properties from the eigenvalues as a tuple,
+        (band gap, cbm, vbm, is_band_gap_direct).
+        """
+
+        vbm = -float("inf")
+        vbm_kpoint = None
+        cbm = float("inf")
+        cbm_kpoint = None
+        for spin, d in self.eigenvalues.items():
+            for k, val in enumerate(d):
+                for (eigenval, occu) in val:
+                    if occu > self.occu_tol and eigenval > vbm:
+                        vbm = eigenval
+                        vbm_kpoint = k
+                    elif occu <= self.occu_tol and eigenval < cbm:
+                        cbm = eigenval
+                        cbm_kpoint = k
+        return max(cbm - vbm, 0), cbm, vbm, vbm_kpoint == cbm_kpoint
+
+
 class Wavederf:
     """
     Object for reading a WAVEDERF file.

pymatgen/io/vasp/tests/test_outputs.py

@@ -20,7 +20,7 @@
 from pymatgen.io.vasp.inputs import Kpoints, Poscar
 from pymatgen.io.vasp.outputs import Chgcar, Locpot, Oszicar, Outcar, \
     Vasprun, Procar, Xdatcar, Dynmat, BSVasprun, UnconvergedVASPWarning, \
-    VaspParserError, Wavecar, Waveder, Elfcar
+    VaspParserError, Wavecar, Waveder, Elfcar, Eigenval
 from pymatgen import Spin, Orbital, Lattice, Structure
 from pymatgen.entries.compatibility import MaterialsProjectCompatibility
 from pymatgen.electronic_structure.core import Magmom
@@ -1412,6 +1412,37 @@ def test_get_parchg(self):
         self.assertFalse(np.all(c.data['total'] > 0.))
 
 
+class EigenvalTest(PymatgenTest):
+    _multiprocess_shared_ = True
+
+    def test_init(self):
+        eig = Eigenval(self.TEST_FILES_DIR / 'EIGENVAL.gz')
+        self.assertEqual(eig.ispin, 1)
+        self.assertEqual(eig.nkpt, len(eig.kpoints))
+        self.assertEqual(eig.nkpt, len(eig.kpoints_weights))
+        self.assertEqual(eig.nkpt, eig.eigenvalues[Spin.up].shape[0])
+        self.assertEqual(eig.nelect, 16)
+        self.assertEqual(eig.nbands, eig.eigenvalues[Spin.up].shape[1])
+        self.assertTrue(np.max(eig.eigenvalues[Spin.up]) > 0)
+        self.assertTrue(np.min(eig.eigenvalues[Spin.up]) < 0)
+
+    def test_ispin2(self):
+        eig = Eigenval(self.TEST_FILES_DIR / 'EIGENVAL.ispin2.gz')
+        self.assertEqual(eig.ispin, 2)
+        self.assertEqual(eig.nkpt, eig.eigenvalues[Spin.up].shape[0])
+        self.assertEqual(eig.nbands, eig.eigenvalues[Spin.up].shape[1])
+        self.assertEqual(eig.nkpt, eig.eigenvalues[Spin.down].shape[0])
+        self.assertEqual(eig.nbands, eig.eigenvalues[Spin.down].shape[1])
+
+    def test_eigenvalue_band_properties(self):
+        eig = Eigenval(self.TEST_FILES_DIR / 'EIGENVAL.gz')
+        props = eig.eigenvalue_band_properties
+        self.assertAlmostEqual(props[0], 6.4153, places=4)
+        self.assertAlmostEqual(props[1], 7.5587, places=4)
+        self.assertAlmostEqual(props[2], 1.1434, places=4)
+        self.assertEqual(props[3], False)
+
+
 class WavederTest(PymatgenTest):
     _multiprocess_shared_ = True