From 2b71180e8e226196187a809a8c6d9c90f30bdca5 Mon Sep 17 00:00:00 2001
From: "Oleg E. Peil" <Oleg.Peil@unige.ch>
Date: Thu, 22 Oct 2015 16:15:49 +0200
Subject: [PATCH] Added calculation of DOS to plotools.py

Added a function that allows one to get the non-interacting projected DOS for
newly generated projectors. The DOS is calculated with analytical tetrahedron
integration added previously.
At the moment, the DOS is generated and output for debugging purposes
after the projectors are generated. Eventually, there should be an
option in the input config file requesting the output of DOS for a given
energy mesh.
---
 python/vasp/plotools.py | 79 +++++++++++++++++++++++++++++++++++++++--
 1 file changed, 76 insertions(+), 3 deletions(-)

diff --git a/python/vasp/plotools.py b/python/vasp/plotools.py
index 8b0e0a89..ef7e0bfe 100644
--- a/python/vasp/plotools.py
+++ b/python/vasp/plotools.py
@@ -1,6 +1,7 @@
 
 import itertools as it
 import numpy as np
+import vasp.atm.c_atm_dos as c_atm_dos
 
 np.set_printoptions(suppress=True)
 
@@ -266,8 +267,8 @@ class ProjectorGroup:
                         p_mat[i1:i2, :nb] = shell.proj_win[ion, isp, ik, :nlm, :nb]
 # Now orthogonalize the obtained block projector
                 p_orth, overl, eig = orthogonalize_projector_matrix(p_mat)
-                print "ik = ", ik
-                print overl.real
+#                print "ik = ", ik
+#                print overl.real
 # Distribute back projectors in the same order
                 for ish in self.ishells:
                     shell = self.shells[ish]
@@ -372,7 +373,7 @@ class ProjectorShell:
 
 ################################################################################
 #
-# select_projectors
+# density_matrix
 #
 ################################################################################
     def density_matrix(self, el_struct, site_diag=True, spin_diag=True):
@@ -405,6 +406,70 @@ class ProjectorShell:
 
         return occ_mats, overlaps
 
+################################################################################
+#
+# density_of_states
+#
+################################################################################
+    def density_of_states(self, el_struct, emesh):
+        """
+        Returns projected DOS for the shell.
+        """
+        nion, ns, nk, nlm, nbtot = self.proj_win.shape
+
+# There is a problem with data storage structure of projectors that will
+# make life more complicated. The problem is that band-indices of projectors
+# for different k-points do not match because we store 'nb_max' values starting
+# from 0.
+        nb_max = self.ib_max - self.ib_min + 1
+        ns_band = self.ib_win.shape[1]
+
+        ne = len(emesh)
+        dos = np.zeros((ne, ns, nion, nlm))
+        w_k = np.zeros((nk, nb_max, ns, nion, nlm), dtype=np.complex128)
+        for isp in xrange(ns):
+            for ik in xrange(nk):
+                is_b = min(isp, ns_band)
+                ib1 = self.ib_win[ik, is_b, 0]
+                ib2 = self.ib_win[ik, is_b, 1] + 1
+                for ib_g in xrange(ib1, ib2):
+                    for io in xrange(nion):
+# Note the difference between 'ib' and 'ibn':
+#  'ib'  counts from 0 to 'nb_k - 1'
+#  'ibn' counts from 'ib1 - ib_min' to 'ib2 - ib_min'
+                        ib = ib_g - ib1
+                        ibn = ib_g - self.ib_min
+                        proj_k = self.proj_win[isp, io, ik, :, ib]
+                        w_k[ik, ib, io, :] = proj_k * proj_k.conj()
+
+#        eigv_ef = el_struct.eigvals[ik, ib, isp] - el_struct.efermi
+        itt = el_struct.kmesh['itet'].T
+# k-indices are starting from 0 in Python
+        itt[1:, :] -= 1
+        for isp in xrange(ns):
+            for ib, eigk in enumerate(el_struct.eigvals[:, self.ib_min:self.ib_max+1, isp].T):
+                for ie, e in enumerate(emesh):
+                    eigk_ef = eigk - el_struct.efermi
+                    cti = c_atm_dos.dos_weights_3d(eigk_ef, e, itt)
+                    for im in xrange(nlm):
+                        for io in xrange(nion):
+                            dos[ie, isp, io, im] += np.sum((cti * w_k[itt[1:, :], ib, isp, io, im].real).sum(0) * itt[0, :])
+
+        dos *= 2 * el_struct.kmesh['volt']
+#        for isp in xrange(ns):
+#            for ik, weight, occ in it.izip(it.count(), kweights, occnums[isp, :, :]):
+#                for io in xrange(nion):
+#                    proj_k = self.proj_win[isp, io, ik, ...]
+#                    occ_mats[isp, io, :, :] += np.dot(proj_k * occ[ib1:ib2],
+#                                                 proj_k.conj().T).real * weight
+#                    overlaps[isp, io, :, :] += np.dot(proj_k,
+#                                                 proj_k.conj().T).real * weight
+
+#        if not symops is None:
+#            occ_mats = symmetrize_matrix_set(occ_mats, symops, ions, perm_map)
+
+        return dos
+
 
 
 ################################################################################
@@ -446,6 +511,14 @@ def generate_plo(conf_pars, el_struct):
         print
         print "Overlap:"
         print ov
+        print
+        print "Evaluating DOS..."
+        emesh = np.linspace(-3.0, 7.0, 4001)
+        dos = pshells[pgroup.ishells[0]].density_of_states(el_struct, emesh)
+        de = emesh[1] - emesh[0]
+        ntot = (dos[1:,...] + dos[:-1,...]).sum(0) / 2 * de
+        print "  Total number of states:", ntot
+        np.savetxt('pdos.dat', np.vstack((emesh.T, dos[:, 0, 0, :].T)).T)
         pgroups.append(pgroup)
 
     return pshells, pgroups