plovasp: documentation and cleaning

python/converters/plovasp/plotools.py

@@ -205,7 +205,7 @@ 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)
+        hk_output(pars, el_struct, pgroups)
 
 
 # TODO: k-points with weights should be stored once and for all
@@ -384,7 +384,6 @@ 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):
@@ -403,47 +402,36 @@ def plo_output(conf_pars, el_struct, pshells, pgroups):
 # plo_output
 #
 ################################################################################
-def hk_output(conf_pars, el_struct, pshells, pgroups):
+def hk_output(conf_pars, el_struct, 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
+    The Hk for each groups is stored in a '<basename>.hk<Ng>' file. The format is
-    defined in the Hk dft_tools format 
+    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
-    ...
 
+    nk                  # number of k-points
+    n_el                # electron density
+    n_sh                # number of total atomic shells
+    at sort l dim       # atom, sort, l, dim
+    at sort l dim       # atom, sort, l, dim
+    n_cr_sh             # number of correlated shells
+    at sort l dim SO irep # atom, sort, l, dim, SO, irep
+    n_irrep dim_irrep   # number of ireps, dim of irep
+    
+    After these header lines, the file has to contain the Hamiltonian matrix
+    in orbital space. The standard convention is that you give for each k-point
+    first the matrix of the real part, then the matrix of the imaginary part,
+    and then move on to the next k-point.
+    
     """
     
-    print 'writing hk'
-    for ig, pgroup in enumerate(pgroups):
-        hk_fname = conf_pars.general['basename'] + '.hk%i'%(ig + 1)
 
+    for ig, pgroup in enumerate(pgroups):
+        
+        hk_fname = conf_pars.general['basename'] + '.hk%i'%(ig + 1)
+        print "  Storing HK-group file '%s'..."%(hk_fname)
 
         head_corr_shells = []
         head_shells = []
@@ -462,32 +450,26 @@ def hk_output(conf_pars, el_struct, pshells, pgroups):
             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
-# Eigenvalues within the window
             nk, nband, ns_band = el_struct.eigvals.shape
-            f.write('%i # number of kpoints\n'%nk)
+            f.write('%i          # number of kpoints\n'%nk)
-            f.write('%f # electron density\n'%1.0)
+            f.write('{0:0.4f}      # electron density\n'.format(pgroup.nelect_window(el_struct)))
-            f.write('%i # number of shells\n'%len(head_shells))
+            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 %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))
+            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('%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')
+            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):

python/converters/plovasp/proj_group.py

@@ -81,6 +81,19 @@ class ProjectorGroup:
         self.ib_min = ib_min
         self.ib_max = ib_max
         self.nb_max = ib_max - ib_min + 1
+        
+        
+        
+        if gr_pars['complement']:
+            n_bands = self.ib_win[:,:,1] - self.ib_win[:,:,0]
+            n_orbs = sum([x.ndim for x in self.shells])
+            assert np.all( n_bands == n_bands[0,0] ), "At each band the same number of bands has to be selected for calculating the complement (to end up with an equal number of orbitals at each k-point)."
+            
+            if n_orbs == n_bands[0,0]+1:
+                gr_pars['complement'] = False
+                print "\nWARNING: The total number of orbitals in this group is  "
+                print "equal to the number of bands. Setting COMPLEMENT to FALSE!\n"
+        
 
 # Select projectors within the energy window
         for ish in self.ishells:
@@ -222,7 +235,7 @@ class ProjectorGroup:
         
         """
 
-        print 'claculating complemet'
+        print '\nCalculating complement\n'
  
         block_maps, ndim = self.get_block_matrix_map()
         _, ns, nk, _, _ = self.shells[0].proj_win.shape
@@ -267,9 +280,6 @@ class ProjectorGroup:
         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
                                                 
                                                 
 ################################################################################

python/converters/plovasp/proj_shell.py

@@ -493,43 +493,14 @@ class ComplementShell(ProjectorShell):
         self.proj_win = proj_compl
         
     def extract_tmatrices(self, sh_pars):
-        pass
+        raise Exception('not implemented')
 
     def local_hamiltonian(self, el_struct, site_diag=True, spin_diag=True):
-        """
+        raise Exception('not implemented')
-        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):
-        """
+        raise Exception('not implemented')
-        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
+        raise Exception('not implemented')
-        ne = len(emesh)
-        dos = np.zeros((ne, ns, nion, nlm))
-        return dos