read eigenvalues, kpoints

parser/parser-siesta/parser.py

@@ -148,6 +148,101 @@ H                     1                    # Species label, number of l-shells
 """
 
 
+def locate_files(dirname, label):
+    files = {}
+    for fileid in ['EIG', 'KP']:
+        path = os.path.join(dirname, '%s.%s' % (label, fileid))
+        if os.path.isfile(path):
+            files[fileid] = path
+        # Warn if files are not present?
+        #
+        # Also: input file
+        #       input parser logfile
+        #
+        # what else?  We already get force/stress/positions from stdout.
+    return files
+
+class SiestaContext(object):
+    def __init__(self):
+        self._is_last_configuration = False
+        self.label = None
+        self.fname = None
+        self.dirname = None  # Base directory of calculations
+        self.files = None
+
+    def set_label(self, parser):
+        # ASSUMPTION: the parser fIn is in the 'root' of whatever was uploaded.
+        # This may not be true.  Figure out how to do this in general.
+        line = parser.fIn.readline()
+        assert line.startswith('reinit: System Label:')
+        self.label = line.split()[-1]
+        self.files = locate_files(self.dirname, self.label)
+
+    def startedParsing(self, fname, parser):
+        self.fname = fname
+        path = os.path.abspath(fname)
+        self.dirname, _ = os.path.split(path)
+
+    #def is_last_configuration(self, parser):
+    #    self._is_last_configuration = True
+
+    def onClose_section_eigenvalues(self, backend, gindex, section):
+        self.read_eigenvalues(backend)
+
+    #def onClose_section_single_configuration_calculation(self, backend, gindex,
+    #                                                     section):
+    #    if self._is_last_configuration:
+    #
+
+    def read_eigenvalues(self, backend):
+        eigfile = self.files.get('EIG')
+        if eigfile is None:
+            return
+
+        with open(eigfile) as fd:
+            eFermi = float(next(fd).strip())
+            nbands, nspins, nkpts = [int(n) for n in next(fd).split()]
+
+            tokens = []
+            for line in fd:
+                tokens.extend(line.split())
+
+        tokens = iter(tokens)
+
+        eps = np.empty((nspins, nkpts, nbands))
+        for k in range(nkpts):
+            kindex = int(next(tokens))
+            assert k + 1 == kindex
+            for s in range(nspins):
+                eps[s, k, :] = [next(tokens) for _ in range(nbands)]
+        unread = list(tokens)
+        assert len(unread) == 0
+        assert s == nspins - 1
+        assert k == nkpts - 1
+
+        # Where does SIESTA store the occupations?
+        backend.addArrayValues('eigenvalues_values', convert_unit(eps, 'eV'))
+
+        kpfile = self.files.get('KP')
+        if kpfile is None:
+            return
+
+
+        with open(kpfile) as fd:
+            tokens = fd.read().split()
+
+        nkpts = int(tokens[0])
+        data = np.array(tokens[1:], object).reshape(-1, 5)
+        print('data', data)
+        coords = data[:, 1:4].astype(float)
+        weights = data[:, 4].astype(float)
+        backend.addArrayValues('eigenvalues_kpoints', coords)
+        # XXX metadata for Fermi level?
+        # k-point weights?
+
+
+context = SiestaContext()
+
 infoFileDescription = SM(
     name='root',
     weak=True,
@@ -160,6 +255,8 @@ infoFileDescription = SM(
     subMatchers=[
         SM(r'Siesta Version: (?P<program_name>siesta)-(?P<program_version>\S+)',
            name='name&version', required=True),
+        SM(r'reinit: System Label:\s*\S*', name='syslabel', forwardMatch=True,
+           adHoc=context.set_label),
         ArraySM(r'siesta: Atomic coordinates \(Bohr\) and species',
                 r'siesta:\s*\S+\s+\S+\s+\S+\s+\S+\s+\S+',
                 add_positions_and_labels),
@@ -216,12 +313,14 @@ infoFileDescription = SM(
                              siesta_energy('Band Struct\.', 'energy_sum_eigenvalues'),
                              siesta_energy('Kinetic', 'electronic_kinetic_energy'),
                              siesta_energy('Hartree', 'energy_electrostatic'),
-                          #siesta_energy('Ext\. field', ''),
+                             #siesta_energy('Ext\. field', ''),
                              siesta_energy('Exch\.-corr\.', 'energy_XC'),
-                          #siesta_energy('Ion-electron', ''),
-                          #siesta_energy('Ion-Ion', ''),
-                          #siesta_energy('Ekinion', ''),
+                             #siesta_energy('Ion-electron', ''),
+                             #siesta_energy('Ion-Ion', ''),
+                             #siesta_energy('Ekinion', ''),
                              siesta_energy('Total', 'energy_total'),
+                             SM(r'', weak=True, name='trigger_readeig',
+                                sections=['section_eigenvalues']),
                          ]),
                       ArraySM(r'siesta: Stress tensor \(static\) \(eV/Ang\*\*3\):',
                               r'siesta:\s*\S+\s+\S+\s+\S+',
@@ -230,18 +329,15 @@ infoFileDescription = SM(
                               r'siesta: Electric dipole \(Debye\)\s*=',
                               get_dipole,
                               forwardMatch=True),
-                   # The purpose of the following matcher is to parse all lines
+                      # The purpose of the following matcher is to parse all lines
                       SM(r'x^', name='end')
                   ])
            ])
     ])
 
-class SiestaContext(object):
-    def startedParsing(self, fname, parser):
-        pass
 
 mainFunction(mainFileDescription=infoFileDescription,
              metaInfoEnv=metaInfoEnv,
              parserInfo=parser_info,
              cachingLevelForMetaName={},
-             superContext=SiestaContext())
+             superContext=context)