From 54b9857aa56de068d5d66e42bbb1649a13c1052d Mon Sep 17 00:00:00 2001
From: "Oleg E. Peil" <Oleg.Peil@unige.ch>
Date: Fri, 16 Oct 2015 16:27:49 +0200
Subject: [PATCH] Added density and overlap matrix output ot ElStruct

The new method in ElectronicStructure allows one to output
denisty and overlap matrices originating from the raw projectors
read from PROJCAR (LOCPROJ). This output is mainly intended for debug purposes.
---
 python/vasp/.gitignore  |  1 +
 python/vasp/elstruct.py | 45 +++++++++++++++++++++++++++++++++++++++++
 python/vasp/main.py     |  1 +
 3 files changed, 47 insertions(+)

diff --git a/python/vasp/.gitignore b/python/vasp/.gitignore
index 0d20b648..442c7f54 100644
--- a/python/vasp/.gitignore
+++ b/python/vasp/.gitignore
@@ -1 +1,2 @@
 *.pyc
+.*.sw*
diff --git a/python/vasp/elstruct.py b/python/vasp/elstruct.py
index d66f9d3e..d7625ac1 100644
--- a/python/vasp/elstruct.py
+++ b/python/vasp/elstruct.py
@@ -74,4 +74,49 @@ class ElectronicStructure:
 # Check that the number of band is the same in PROJCAR and EIGENVAL
         assert nb_plo == self.nband, "PLOCAR is inconsistent with EIGENVAL (number of bands)"
 
+    def debug_density_matrix(self):
+        """
+        Calculate and output the density and overlap matrix out of projectors defined in el_struct.
+        """
+        plo = self.proj_raw
+        nproj, ns, nk, nb = plo.shape
+        ions = list(set([param['isite'] for param in self.proj_params]))
+        nions = len(ions)
+        norb = nproj / nions
+
+        den_mat = np.zeros((ns, nproj, nproj), dtype=np.float64)
+        overlap = np.zeros((ns, nproj, nproj), dtype=np.float64)
+#        ov_min = np.ones((ns, nproj, nproj), dtype=np.float64) * 100.0
+#        ov_max = np.zeros((ns, nproj, nproj), dtype=np.float64)
+        for ispin in xrange(ns):
+            for ik in xrange(nk):
+                kweight = self.kmesh['kweights'][ik]
+                occ = self.ferw[ispin, ik, :]
+                den_mat[ispin, :, :] += np.dot(plo[:, ispin, ik, :] * occ, plo[:, ispin, ik, :].T.conj()).real * kweight
+                ov = np.dot(plo[:, ispin, ik, :], plo[:, ispin, ik, :].T.conj()).real
+                overlap[ispin, :, :] += ov * kweight
+#                ov_max = np.maximum(ov, ov_max)
+#                ov_min = np.minimum(ov, ov_min)
+
+# Output only the site-diagonal parts of the matrices
+        for ispin in xrange(ns):
+            print
+            print "  Spin:", ispin + 1
+            for io, ion in enumerate(ions):
+                print "  Site:", ion
+                iorb_inds = [(ip, param['m']) for ip, param in enumerate(self.proj_params) if param['isite'] == ion]
+                norb = len(iorb_inds)
+                dm = np.zeros((norb, norb))
+                ov = np.zeros((norb, norb))
+                for ind, iorb in iorb_inds:
+                    dm[iorb, :] = den_mat[ispin, ind, :]
+                    ov[iorb, :] = overlap[ispin, ind, :]
+
+                print "  Density matrix" + (12*norb - 12)*" " + "Overlap"
+                for drow, dov in zip(dm, ov):
+                    out = ''.join(map("{0:12.7f}".format, drow))
+                    out += "    "
+                    out += ''.join(map("{0:12.7f}".format, dov))
+                    print out
+
 
diff --git a/python/vasp/main.py b/python/vasp/main.py
index 665fdbc3..4d1fa933 100644
--- a/python/vasp/main.py
+++ b/python/vasp/main.py
@@ -23,5 +23,6 @@ if __name__ == '__main__':
     pars.parse_input()
     vasp_data = vaspio.VaspData(vasp_dir)
     el_struct = ElectronicStructure(vasp_data)
+    el_struct.debug_density_matrix()
     pshells, pgroups = generate_plo(pars, el_struct)
     output_as_text(pars, el_struct, pshells, pgroups)