From 0eb574b5c8e28c0bec30da462c166ff10d13b775 Mon Sep 17 00:00:00 2001
From: "Oleg E. Peil" <Oleg.Peil@unige.ch>
Date: Tue, 10 Nov 2015 16:40:46 +0100
Subject: [PATCH] Added DOSMESH option to section [General]

If option DOSMESH is specified a projected DOS for each shell
will be output. Energy mesh parameters are given in DOSMESH as

DOSMESH = [EMIN  EMAX]  N_POINTS

The parameters in the brackets [] are optional. If only the number
of points is specified the energy range is taken to be the same
as the projection energy window.
---
 python/vasp/inpconf.py  | 44 ++++++++++++++++++++++++++++++++++++++++-
 python/vasp/plotools.py | 28 ++++++++++++++++++--------
 2 files changed, 63 insertions(+), 9 deletions(-)

diff --git a/python/vasp/inpconf.py b/python/vasp/inpconf.py
index e27d082c..e0b8fafc 100644
--- a/python/vasp/inpconf.py
+++ b/python/vasp/inpconf.py
@@ -70,7 +70,8 @@ class ConfigParameters:
 
         self.gen_optional = {
             'basename' : ('basename', str, 'vasp'),
-            'efermi' : ('efermi', float)}
+            'efermi' : ('efermi', float),
+            'dosmesh': ('dosmesh', self.parse_string_dosmesh)}
 
 #
 # Special parsers
@@ -97,6 +98,7 @@ class ConfigParameters:
             i1, i2 = tuple(map(int, match.groups()[:2]))
             mess = "First index of the range must be smaller or equal to the second"
             assert i1 <= i2, mess
+# Note that we need to subtract 1 from VASP indices
             ion_list = np.array(range(i1 - 1, i2))
         else:
 # Check if a set of indices is given
@@ -165,6 +167,46 @@ class ConfigParameters:
 
         return mat
 
+################################################################################
+#
+# parse_string_ion_list()
+#
+################################################################################
+    def parse_string_dosmesh(self, par_str):
+        """
+        Two formats are accepted:
+
+          # Two floats (energy range) and an integer (number of energy points).
+
+          # One integer (number of energy points). In this case the energy
+            range is taken to be equal to EMIN, EMAX of a shell.
+
+        The parser returns a dictionary:
+          {'n_points': int,
+           'emin': float,
+           'emax': float}
+
+        If the second option is used, 'emin' and 'emax' are undefined
+        and set to 'nan'.
+        """
+        stmp = par_str.split()
+        if len(stmp) == 3:
+            emin, emax = float(stmp[0]), float(stmp[1])
+            n_points = int(stmp[2])
+        elif len(stmp) == 1:
+            n_points = int(stmp[0])
+            emin = emax = float('nan')
+        else:
+            err_mess = "DOSMESH must be either 'EMIN EMAX NPOINTS' or 'NPOINTS'"
+            raise ValueError(err_mess)
+
+        dos_pars = {
+            'n_points': n_points,
+            'emin': emin,
+            'emax': emax}
+
+        return dos_pars
+
 ################################################################################
 #
 # parse_parameter_set()
diff --git a/python/vasp/plotools.py b/python/vasp/plotools.py
index f0ad76cd..f7e11db0 100644
--- a/python/vasp/plotools.py
+++ b/python/vasp/plotools.py
@@ -520,14 +520,26 @@ 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)
+        if 'dosmesh' in conf_pars.general:
+            print
+            print "Evaluating DOS..."
+            mesh_pars = conf_pars.general['dosmesh']
+            if np.isnan(mesh_pars['emin']):
+                dos_emin = pgroup.emin
+                dos_emax = pgroup.emax
+            else:
+                dos_emin = mesh_pars['emin']
+                dos_emax = mesh_pars['emax']
+            n_points = mesh_pars['n_points']
+
+            emesh = np.linspace(dos_emin, dos_emax, n_points)
+            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
+            for io in xrange(dos.shape[2]):
+                np.savetxt('pdos_%i.dat'%(io), np.vstack((emesh.T, dos[:, 0, io, :].T)).T)
+
         pgroups.append(pgroup)
 
     return pshells, pgroups