From 66fac2f1bd8288cb19c53cad38783a53f9d289b2 Mon Sep 17 00:00:00 2001
From: "Oleg E. Peil" <Oleg.Peil@unige.ch>
Date: Wed, 14 Oct 2015 15:58:45 +0200
Subject: [PATCH] Added preliminary PROJCAR parser to vaspio

This python-parser is a prototype of a future parser that will probably
be using only LOCPROJ (which is going to be modified).
At the moment, one has to use the first line of LOCPROJ to determine
the array dimensions and parse PROJCAR because it contains relevant information
on projectors (such as site and orbital character).

Note that in the previous implementation relying on the binary PLOCAR-file
the Fermi weights were taken from PLOCAR. In the current version of VASP
(>=5.4.1) the Fermi weights can read in from EIGENVAL.
---
 python/vasp/vaspio.py | 94 ++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 92 insertions(+), 2 deletions(-)

diff --git a/python/vasp/vaspio.py b/python/vasp/vaspio.py
index 1f7ab77c..404c603f 100644
--- a/python/vasp/vaspio.py
+++ b/python/vasp/vaspio.py
@@ -1,6 +1,7 @@
 
 import numpy as np
-import plocar_io.c_plocar_io as c_plocar_io
+import re
+#import plocar_io.c_plocar_io as c_plocar_io
 
 def read_lines(filename):
     r"""
@@ -77,7 +78,96 @@ class Plocar:
         if vasp_dir[-1] != '/':
             vasp_dir += '/'
 
-        self.params, self.plo, self.ferw = c_plocar_io.read_plocar(vasp_dir + plocar_filename)
+#        self.params, self.plo, self.ferw = c_plocar_io.read_plocar(vasp_dir + plocar_filename)
+        self.params, self.plo = self.temp_parser(projcar_filename=vasp_dir + "PROJCAR", locproj_filename=vasp_dir + "LOCPROJ")
+
+    def temp_parser(self, projcar_filename='PROJCAR', locproj_filename='LOCPROJ'):
+        r"""
+        Parses PROJCAR (and partially LOCPROJ) to get VASP projectors.
+
+        This is a prototype parser that should eventually be written in C for
+        better performance on large files.
+
+        Returns projector parameters (site/orbital indices etc.) and an array
+        with projectors.
+        """
+        orb_labels = ["s", "pz", "px", "py", "dz2", "dxz", "dyz", "dx2-y2", "dxy",
+                      "fz3", "fxz2", "fyz2", "fz(x2-y2)", "fxyz", "fx(x2-3y2)", "fy(3x2-y2)"]
+
+        def lm_to_l_m(lm):
+            l = int(np.sqrt(lm))
+            m = lm - l*l
+            return l, m
+
+# Read the first line of LOCPROJ to get the dimensions
+        with open(locproj_filename, 'rt') as f:
+            line = f.readline()
+            nproj, nspin, nk, nband = map(int, line.split())
+
+        plo = np.zeros((nproj, nspin, nk, nband), dtype=np.complex128)
+        params = [{} for i in xrange(nproj)]
+
+        iproj_site = 0
+        is_first_read = True
+        with open(projcar_filename, 'rt') as f:
+            line = self.search_for(f, "^ *ISITE")
+            while line:
+                isite = int(line.split()[1])
+                if not is_first_read:
+                    for il in xrange(norb):
+                        params[iproj_site * norb + il]['isite'] = isite
+                        params[iproj_site * norb + il]['l'] = params[(iproj_site - 1) * norb + il]['l']
+                        params[iproj_site * norb + il]['m'] = params[(iproj_site - 1) * norb + il]['m']
+
+                for ispin in xrange(nspin):
+                    for ik in xrange(nk):
+# Parse the orbital labels and convert them to l,m-indices
+                        line = self.search_for(f, "^ *band")
+                        if is_first_read:
+                            cpatt = re.compile("lm= *([^\s]+)")
+                            labels = re.findall(cpatt, line)
+                            norb = len(labels)
+                            for il, label in enumerate(labels):
+                                lm = orb_labels.index(label)
+                                l, m = lm_to_l_m(lm)
+
+# For the first read 'iproj_site = 0' and only orbital index 'il' is used
+                                params[il]['isite'] = isite
+                                params[il]['l'] = l
+                                params[il]['m'] = m
+
+                            is_first_read = False
+
+# Read the block of nk * ns * nband complex numbers
+                        for ib in xrange(nband):
+                            line = f.readline()
+                            rtmp = map(float, line.split()[1:])
+                            for il in xrange(norb):
+                                ctmp = complex(rtmp[2 * il], rtmp[2 * il + 1])
+                                plo[iproj_site * norb + il, ispin, ik, ib] = ctmp
+
+# End of site-block
+                iproj_site += 1
+                line = self.search_for(f, "^ *ISITE")
+
+        print "Read parameters:"
+        for il, par in enumerate(params):
+            print il, " -> ", par
+
+        return params, plo
+
+
+    def search_for(self, f, patt):
+        r"""
+        Reads file 'f' until pattern 'patt' is encountered and returns
+        the corresponding line.
+        """
+        cpatt = re.compile(patt)
+        line = "x"
+        while not re.match(cpatt, line) and line:
+            line = f.readline()
+
+        return line
 
 
 ################################################################################