diff --git a/python/converters/plovasp/plotools.py b/python/converters/plovasp/plotools.py
index e980858d..000c2123 100644
--- a/python/converters/plovasp/plotools.py
+++ b/python/converters/plovasp/plotools.py
@@ -38,6 +38,7 @@ import itertools as it
 import numpy as np
 from proj_group import ProjectorGroup
 from proj_shell import ProjectorShell
+from proj_shell import ComplementShell
 
 np.set_printoptions(suppress=True)
 
@@ -124,24 +125,29 @@ def generate_plo(conf_pars, el_struct):
     for gr_par in conf_pars.groups:
         pgroup = ProjectorGroup(gr_par, pshells, eigvals)
         pgroup.orthogonalize()
+        if gr_par['complement']:
+            pgroup.complement(eigvals)
+        if conf_pars.general['hk']:
+            pgroup.calc_hk(eigvals)
 # DEBUG output
         print "Density matrix:"
         nimp = 0.0
         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
+            if  not isinstance(pshells[pgroup.ishells[ish]],ComplementShell):
+                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
@@ -152,12 +158,13 @@ def generate_plo(conf_pars, el_struct):
         print
         print "Local Hamiltonian:"
         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))
+            if  not isinstance(pshells[pgroup.ishells[ish]],ComplementShell):
+                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
@@ -173,13 +180,14 @@ def generate_plo(conf_pars, el_struct):
 
             emesh = np.linspace(dos_emin, dos_emax, n_points)
             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)
+                if  not isinstance(pshells[pgroup.ishells[ish]],ComplementShell):
+                    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)
 
@@ -196,6 +204,8 @@ def output_as_text(pars, el_struct, pshells, pgroups):
     """
     ctrl_output(pars, el_struct, len(pgroups))
     plo_output(pars, el_struct, pshells, pgroups)
+    if pars.general['hk']:
+        hk_output(pars, el_struct, pshells, pgroups)
 
 
 # TODO: k-points with weights should be stored once and for all
@@ -388,4 +398,102 @@ def plo_output(conf_pars, el_struct, pshells, pgroups):
                                     f.write("{0:16.10f}{1:16.10f}\n".format(p.real, p.imag))
                                 f.write("\n")
 
+################################################################################
+#
+# plo_output
+#
+################################################################################
+def hk_output(conf_pars, el_struct, pshells, pgroups):
+    """
+    Outputs HK into text file.
 
+    Filename is defined by <basename> that is passed from config-file.
+
+    The Hk for all groups is stored in a '<basename>.hk' file. The format is as
+    defined in the Hk dft_tools format 
+
+    # Energy window: emin, emax
+    ib_min, ib_max
+    nelect
+    # Eigenvalues
+    isp, ik1, kx, ky, kz, kweight
+    ib1, ib2
+    eig1
+    eig2
+    ...
+    eigN
+    ik2, kx, ky, kz, kweight
+    ...
+
+    # Projected shells
+    Nshells
+    # Shells: <shell indices>
+    # Shell <1>
+    Shell 1
+    ndim
+    # complex arrays: plo(ns, nion, ndim, nb)
+    ...
+    # Shells: <shell indices>
+    # Shell <2>
+    Shell 2
+    ...
+
+    """
+    
+    print 'writing hk'
+    for ig, pgroup in enumerate(pgroups):
+        hk_fname = conf_pars.general['basename'] + '.hk%i'%(ig + 1)
+
+
+        head_corr_shells = []
+        head_shells = []
+        for ish in pgroup.ishells:
+
+            shell = pgroup.shells[ish]
+            
+            sh_dict = {}
+            sh_dict['shell_index'] = ish
+            sh_dict['lorb'] = shell.lorb
+            sh_dict['ndim'] = shell.ndim
+# 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]
+# Derive sorts from equivalence classes
+            sh_dict['ion_list'] = ion_output
+            sh_dict['ion_sort'] = shell.ion_sort
+
+# TODO: add the output of transformation matrices
+
+            head_shells.append(sh_dict)
+            if shell.corr:
+                head_corr_shells.append(sh_dict)
+
+        #head_dict['shells'] = head_shells
+
+        #header = json.dumps(head_dict, indent=4, separators=(',', ': '))
+
+        with open(hk_fname, 'wt') as f:
+                       
+# Eigenvalues within the window
+            nk, nband, ns_band = el_struct.eigvals.shape
+            f.write('%i # number of kpoints\n'%nk)
+            f.write('%f # electron density\n'%1.0)
+            f.write('%i # number of shells\n'%len(head_shells))
+            for head in head_shells:
+                f.write('%i %i %i %i # atom sort l dim\n'%(head['ion_list'][0],head['ion_sort'][0],head['lorb'],head['ndim']))
+            f.write('%i # number of correlated shells\n'%len(head_corr_shells))
+            for head in head_corr_shells:
+                f.write('%i %i %i %i 0 0 # atom sort l dim SO irrep\n'%(head['ion_list'][0],head['ion_sort'][0],head['lorb'],head['ndim']))
+            f.write('1 5 # number of ireps, dim of irep\n')
+            norbs = pgroup.hk.shape[2]
+            for isp in xrange(ns_band):
+                #f.write("# is = %i\n"%(isp + 1))
+                for ik in xrange(nk):
+                    for io in xrange(norbs):
+                        for iop in xrange(norbs):
+                            f.write(" {0:14.10f}".format(pgroup.hk[isp,ik,io,iop].real))
+                        f.write("\n")
+                    for io in xrange(norbs):
+                        for iop in xrange(norbs):
+                            f.write(" {0:14.10f}".format(pgroup.hk[isp,ik,io,iop].imag))
+                        f.write("\n")
diff --git a/python/converters/plovasp/proj_group.py b/python/converters/plovasp/proj_group.py
index 1ad99dd6..1fff7a04 100644
--- a/python/converters/plovasp/proj_group.py
+++ b/python/converters/plovasp/proj_group.py
@@ -30,7 +30,7 @@ r"""
     Storage and manipulation of projector groups.
 """
 import numpy as np
-
+from proj_shell import ComplementShell
 np.set_printoptions(suppress=True)
 
 ################################################################################
@@ -165,7 +165,113 @@ class ProjectorGroup:
                         nlm = i2 - i1 + 1
                         shell = self.shells[ish]
                         shell.proj_win[ion, isp, ik, :nlm, :nb] = p_orth[i1:i2, :nb]
+                        
+                        
+################################################################################
+#
+# calc_hk
+#
+################################################################################
+    def calc_hk(self, eigvals):
+        """
+        Calculate H(k) for a group by applying the projectors to the eigenvalues.
+        
+        """
 
+        block_maps, ndim = self.get_block_matrix_map()
+
+        _, ns, nk, _, _ = self.shells[0].proj_win.shape
+        p_mat = np.zeros((ndim, self.nb_max), dtype=np.complex128)
+        #print(p_mat.shape)
+        self.hk = np.zeros((ns,nk,ndim,ndim), dtype=np.complex128)
+# Note that 'ns' and 'nk' are the same for all shells
+        for isp in xrange(ns):
+            for ik in xrange(nk):
+                bmin = self.ib_win[ik, isp, 0]
+                bmax = self.ib_win[ik, isp, 1]+1
+                
+                nb = bmax - bmin 
+                #print(bmin,bmax,nb)
+# Combine all projectors of the group to one block projector
+                for bl_map in block_maps:
+                    p_mat[:, :] = 0.0j  # !!! Clean-up from the last k-point and block!
+                    for ibl, block in enumerate(bl_map):
+                        i1, i2 = block['bmat_range']
+                        ish, ion = block['shell_ion']
+                        nlm = i2 - i1 + 1
+                        shell = self.shells[ish]
+                        p_mat[i1:i2, :nb] = shell.proj_win[ion, isp, ik, :nlm, :nb]
+                    #print(p_mat.shape, eigvals[ik,:,isp].shape)
+                    self.hk[isp,ik,:,:] = np.einsum('ij,jk->ik',p_mat*eigvals[ik,bmin:bmax,isp],
+                                                p_mat.transpose().conjugate())
+                    #print(np.linalg.eigvals(self.hk[isp,ik,...]))
+                    #print(eigvals[ik,bmin:bmax,isp])
+                    #print(self.hk.shape)
+                    #self.hk[isp,ik,:,:] = np.dot(p_mat.conjugate()*eigvals[ik,bmin:bmax,isp],
+                    #                            p_mat.transpose())
+
+
+################################################################################
+#
+# complement
+#
+################################################################################
+    def complement(self,eigvals):
+        """
+        Calculate the complement for a group of projectors.
+        
+        """
+
+        print 'claculating complemet'
+ 
+        block_maps, ndim = self.get_block_matrix_map()
+        _, ns, nk, _, _ = self.shells[0].proj_win.shape
+        p_mat = np.zeros((ndim, self.nb_max), dtype=np.complex128)
+        p_full = np.zeros((1,ns,nk,self.nb_max, self.nb_max), dtype=np.complex128)        
+        
+        #print(p_mat.shape)
+        self.hk = np.zeros((ns,nk,ndim,ndim), dtype=np.complex128)
+# Note that 'ns' and 'nk' are the same for all shells
+        orbs_done = 1*ndim
+        while orbs_done < self.nb_max:
+            
+            for isp in xrange(ns):
+                for ik in xrange(nk):
+                    bmin = self.ib_win[ik, isp, 0]
+                    bmax = self.ib_win[ik, isp, 1]+1
+                    
+                    nb = bmax - bmin 
+                    #print(bmin,bmax,nb)
+    # Combine all projectors of the group to one block projector
+                    for bl_map in block_maps:
+                        p_mat[:, :] = 0.0j  # !!! Clean-up from the last k-point and block!
+                        for ibl, block in enumerate(bl_map):
+                            i1, i2 = block['bmat_range']
+                            ish, ion = block['shell_ion']
+                            nlm = i2 - i1 + 1
+                            shell = self.shells[ish]
+                            p_mat[i1:i2, :nb] = shell.proj_win[ion, isp, ik, :nlm, :nb]
+                    p_full[0,isp,ik,:ndim,:] = p_mat
+                    proj_work = p_full[0,isp,ik,:,:] 
+                    overlap = np.dot(proj_work.transpose().conjugate(),proj_work)
+                    work = np.eye(self.nb_max) - overlap
+                    norm = np.sqrt(np.sum(work*work.transpose().conjugate(),axis=1))
+                    max_ind = np.argmax(norm)
+                    p_full[0,isp,ik,orbs_done,:] = work[max_ind]/norm[max_ind]
+                    
+            orbs_done += 1
+        sh_pars = {}
+        sh_pars['lshell'] = -1
+        sh_pars['ions'] = {'nion':1,'ion_list':[[1]]}
+        sh_pars['user_index'] = 'complement'
+        sh_pars['corr']  = False
+        self.shells.append(ComplementShell(sh_pars,p_full[:,:,:,ndim:,:],False))
+        self.ishells.append(self.ishells[-1]+1)
+        print(self.shells)
+        print(self.ishells)
+        #p_full -= p_mat
+                                                
+                                                
 ################################################################################
 #
 # gen_block_matrix_map
diff --git a/python/converters/plovasp/proj_shell.py b/python/converters/plovasp/proj_shell.py
index 4aa6d8ed..463544a0 100644
--- a/python/converters/plovasp/proj_shell.py
+++ b/python/converters/plovasp/proj_shell.py
@@ -450,5 +450,86 @@ class ProjectorShell:
 
         return dos
 
+################################################################################
+################################################################################
+#
+# class ProjectorShell
+#
+################################################################################
+################################################################################
+class ComplementShell(ProjectorShell):
+    """
+    Container of projectors related to a complement shell.
 
 
+    Parameters:
+
+    - sh_pars (dict) : shell parameters from the config-file
+    - proj_compl (numpy.array) : array of complement projectors
+
+    """
+    def __init__(self, sh_pars, proj_compl, nc_flag):
+        self.lorb = sh_pars['lshell']
+        self.ions = sh_pars['ions']
+        self.user_index = sh_pars['user_index']
+        self.corr = sh_pars['corr']
+        self.nc_flag = nc_flag        
+
+
+        #self.lm1 = self.lorb**2
+        #self.lm2 = (self.lorb+1)**2
+
+        self.nion = self.ions['nion']
+# Extract ion list and equivalence classes (ion sorts)
+        self.ion_list = sorted(it.chain(*self.ions['ion_list']))
+        self.ion_sort = []
+        for ion in self.ion_list:
+            for icl, eq_cl in enumerate(self.ions['ion_list']):
+                if ion in eq_cl:
+                    self.ion_sort.append(icl + 1) # Enumerate classes starting from 1
+                    break
+
+        self.ndim = proj_compl.shape[3]
+        self.proj_win = proj_compl
+        
+    def extract_tmatrices(self, sh_pars):
+        pass
+
+    def local_hamiltonian(self, el_struct, site_diag=True, spin_diag=True):
+        """
+        Returns occupation matrix/matrices for the shell.
+        """
+        nion, ns, nk, nlm, nbtot = self.proj_win.shape
+
+        assert site_diag, "site_diag = False is not implemented"
+        assert spin_diag, "spin_diag = False is not implemented"
+
+        loc_ham = np.zeros((ns, nion, nlm, nlm), dtype=np.float64)
+
+        return loc_ham
+        
+    def density_matrix(self, el_struct, site_diag=True, spin_diag=True):
+        """
+        Returns occupation matrix/matrices for the shell.
+        """
+        nion, ns, nk, nlm, nbtot = self.proj_win.shape
+
+#        assert site_diag, "site_diag = False is not implemented"
+        assert spin_diag, "spin_diag = False is not implemented"
+
+        if site_diag:
+            occ_mats = np.zeros((ns, nion, nlm, nlm), dtype=np.float64)
+            overlaps = np.zeros((ns, nion, nlm, nlm), dtype=np.float64)
+        else:
+            ndim = nion * nlm
+            occ_mats = np.zeros((ns, 1, ndim, ndim), dtype=np.float64)
+            overlaps = np.zeros((ns, 1, ndim, ndim), dtype=np.float64)
+
+        return occ_mats, overlaps
+
+    def density_of_states(self, el_struct, emesh):
+        nion, ns, nk, nlm, nbtot = self.proj_win.shape
+        ne = len(emesh)
+        dos = np.zeros((ne, ns, nion, nlm))
+        return dos
+    
\ No newline at end of file