diff --git a/python/converters/plovasp/inpconf.py b/python/converters/plovasp/inpconf.py
index aad926de..9feac451 100644
--- a/python/converters/plovasp/inpconf.py
+++ b/python/converters/plovasp/inpconf.py
@@ -84,7 +84,8 @@ class ConfigParameters:
 
         self.sh_optional = {
             'transform': ('tmatrix', lambda s: self.parse_string_tmatrix(s, real=True)),
-            'transfile': ('tmatrices', self.parse_file_tmatrix)}
+            'transfile': ('tmatrices', self.parse_file_tmatrix),
+            'corr': ('corr', self.parse_string_logical, True)}
 
         self.gr_required = {
             'shells': ('shells', lambda s: map(int, s.split())),
@@ -92,12 +93,16 @@ class ConfigParameters:
 
         self.gr_optional = {
             'normalize' : ('normalize', self.parse_string_logical, True),
-            'normion' : ('normion', self.parse_string_logical, True)}
+            'normion' : ('normion', self.parse_string_logical, True),
+            'complement' : ('complement', self.parse_string_logical, False),
+            'bands': ('bands', self.parse_band_window)}
+
 
         self.gen_optional = {
             'basename' : ('basename', str, 'vasp'),
             'efermi' : ('efermi', float),
-            'dosmesh': ('dosmesh', self.parse_string_dosmesh)}
+            'dosmesh': ('dosmesh', self.parse_string_dosmesh),
+            'hk': ('hk', self.parse_string_logical, False)}
 
 #
 # Special parsers
@@ -205,6 +210,21 @@ class ConfigParameters:
         assert ftmp[0] < ftmp[1], "The first float in EWINDOW must be smaller than the second one"
         return tuple(ftmp)
 
+################################################################################
+#
+# parse_band_window()
+#
+################################################################################
+    def parse_band_window(self, par_str):
+        """
+        Band window is given by two ints, with the first one being smaller
+        than the second one.
+        """
+        ftmp = map(int, par_str.split())
+        assert len(ftmp) == 2, "BWINDOW must be specified by exactly two ints"
+        assert ftmp[0] < ftmp[1], "The first int in BWINDOW must be smaller than the second one"
+        return tuple(ftmp)
+
 ################################################################################
 #
 # parse_string_tmatrix()
diff --git a/python/converters/plovasp/plotools.py b/python/converters/plovasp/plotools.py
index e698d49f..e980858d 100644
--- a/python/converters/plovasp/plotools.py
+++ b/python/converters/plovasp/plotools.py
@@ -117,6 +117,7 @@ def generate_plo(conf_pars, el_struct):
         print "    Orbital l     : %i"%(pshell.lorb)
         print "    Number of ions: %i"%(pshell.nion)
         print "    Dimension     : %i"%(pshell.ndim)
+        print "    Correlated    : %r"%(pshell.corr)
         pshells.append(pshell)
 
     pgroups = []
@@ -125,31 +126,38 @@ def generate_plo(conf_pars, el_struct):
         pgroup.orthogonalize()
 # DEBUG output
         print "Density matrix:"
-        dm_all, ov_all = pshells[pgroup.ishells[0]].density_matrix(el_struct)
         nimp = 0.0
-        spin_fac = 2 if dm_all.shape[0] == 1 else 1
-        for io in xrange(dm_all.shape[1]):
-            print "  Site %i"%(io + 1)
-            dm = spin_fac * dm_all[:, io, : ,:].sum(0)
-            for row in dm:
-                print ''.join(map("{0:12.7f}".format, row))
-            ndm = dm.trace()
-            nimp += ndm
-            print "    trace: ", ndm
+        ov_all = []
+        for ish in pgroup.ishells:
+            print "  Shell %i"%(ish + 1)
+            dm_all, ov_all_ = pshells[ish].density_matrix(el_struct)
+            ov_all.append(ov_all_[0])
+            spin_fac = 2 if dm_all.shape[0] == 1 else 1
+            for io in xrange(dm_all.shape[1]):
+                print "    Site %i"%(io + 1)
+                dm = spin_fac * dm_all[:, io, : ,:].sum(0)
+                for row in dm:
+                    print ''.join(map("{0:14.7f}".format, row))
+                ndm = dm.trace()
+                if pshells[ish].corr:
+                    nimp += ndm
+                print "      trace: ", ndm
         print
         print "  Impurity density:", nimp
         print
         print "Overlap:"
-        for io, ov in enumerate(ov_all[0]):
+        for io, ov in enumerate(ov_all):
             print "  Site %i"%(io + 1)
-            print ov
+            print ov[0,...]
         print
         print "Local Hamiltonian:"
-        loc_ham = pshells[pgroup.ishells[0]].local_hamiltonian(el_struct)
-        for io in xrange(loc_ham.shape[1]):
-            print "  Site %i"%(io + 1)
-            for row in loc_ham[:, io, :, :].sum(0):
-                print ''.join(map("{0:12.7f}".format, row))
+        for ish in pgroup.ishells:
+            print "  Shell %i"%(ish + 1)
+            loc_ham = pshells[pgroup.ishells[ish]].local_hamiltonian(el_struct)
+            for io in xrange(loc_ham.shape[1]):
+                print "    Site %i"%(io + 1)
+                for row in loc_ham[:, io, :, :].sum(0):
+                    print ''.join(map("{0:14.7f}".format, row))
 # END DEBUG output
         if 'dosmesh' in conf_pars.general:
             print
@@ -164,12 +172,14 @@ def generate_plo(conf_pars, el_struct):
             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)
+            for ish in pgroup.ishells:
+                print "  Shell %i"%(ish + 1)
+                dos = pshells[pgroup.ishells[ish]].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_%i.dat'%(ish,io), np.vstack((emesh.T, dos[:, 0, io, :].T)).T)
 
         pgroups.append(pgroup)
 
@@ -322,6 +332,7 @@ def plo_output(conf_pars, el_struct, pshells, pgroups):
             sh_dict['shell_index'] = ish
             sh_dict['lorb'] = shell.lorb
             sh_dict['ndim'] = shell.ndim
+            sh_dict['corr'] = shell.corr
 # Convert ion indices from the internal representation (starting from 0)
 # to conventional VASP representation (starting from 1)
             ion_output = [io + 1 for io in shell.ion_list]
@@ -363,6 +374,7 @@ def plo_output(conf_pars, el_struct, pshells, pgroups):
                 f.write("# Shell %i\n"%(ish))
 
                 nion, ns, nk, nlm, nb = shell.proj_win.shape
+                print(nlm)
                 for isp in xrange(ns):
                     f.write("# is = %i\n"%(isp + 1))
                     for ik in xrange(nk):
diff --git a/python/converters/plovasp/proj_group.py b/python/converters/plovasp/proj_group.py
index ddd49e8d..1ad99dd6 100644
--- a/python/converters/plovasp/proj_group.py
+++ b/python/converters/plovasp/proj_group.py
@@ -67,6 +67,16 @@ class ProjectorGroup:
 
 # Determine the minimum and maximum band numbers
         ib_win, ib_min, ib_max = self.select_bands(eigvals)
+        if 'bands' in gr_pars:
+            nk, nband, ns_band = eigvals.shape
+            ib_win = np.zeros((nk, ns_band, 2), dtype=np.int32)
+            ib_win[:,:,0] = gr_pars['bands'][0]
+            ib_win[:,:,1] = gr_pars['bands'][1]
+            ib_min = gr_pars['bands'][0]
+            ib_max = gr_pars['bands'][1]
+            
+        else:
+            ib_win, ib_min, ib_max = self.select_bands(eigvals)
         self.ib_win = ib_win
         self.ib_min = ib_min
         self.ib_max = ib_max
diff --git a/python/converters/plovasp/proj_shell.py b/python/converters/plovasp/proj_shell.py
index 93c7ae0f..4aa6d8ed 100644
--- a/python/converters/plovasp/proj_shell.py
+++ b/python/converters/plovasp/proj_shell.py
@@ -72,7 +72,8 @@ class ProjectorShell:
         self.lorb = sh_pars['lshell']
         self.ions = sh_pars['ions']
         self.user_index = sh_pars['user_index']
-        self.nc_flag = nc_flag
+        self.corr = sh_pars['corr']
+        self.nc_flag = nc_flag        
 #        try:
 #            self.tmatrix = sh_pars['tmatrix']
 #        except KeyError: