Changed order of indices of 'proj_arr' array

When a ProjectorShell is created it creates a view of the full
projector array with orbital 'ilm' and band 'ib' indices interchanged.
The reason for this is that this corresponds more naturally to the
definition of the projector P_{m\nu} and also allows for multiplications
of projector matrices without additional transposition.

The tests have been modified accordingly.

python/converters/vasp/python/plotools.py

@@ -202,7 +202,7 @@ class ProjectorGroup:
                 i1_bl = i2_bl
             bl_map[ish]['bmat_blocks'] = bmat_bl
 
-        ndim = i2
+        ndim = i2_bl
         p_mat = np.zeros((ndim, nb_max), dtype=np.complex128)
         for isp in xrange(ns):
             for ik in xrange(nk):
@@ -257,7 +257,9 @@ class ProjectorShell:
         self.lm2 = (self.lorb+1)**2
 
 # Pre-select a subset of projectors (this should be an array view => no memory is wasted)
-        self.proj_arr = proj_raw[self.ion_list, :, :, :, self.lm1:self.lm2]
+# !!! This sucks but I have to change the order of 'ib' and 'ilm' indices here
+# This should perhaps be done right after the projector array is read from PLOCAR
+        self.proj_arr = proj_raw[self.ion_list, :, :, :, self.lm1:self.lm2].transpose((0, 1, 2, 4, 3))
 
 ################################################################################
 #
@@ -274,9 +276,9 @@ class ProjectorShell:
 
 # Set the dimensions of the array
         nb_win = self.nb_max - self.nb_min + 1
-        nion, ns, nk, nbtot, nlm = self.proj_arr.shape
+        nion, ns, nk, nlm, nbtot = self.proj_arr.shape
 # !!! Note that the order is changed below !!!
-        self.proj_win = np.zeros((nion, ns, nk, nb_win, nlm), dtype=np.complex128)
+        self.proj_win = np.zeros((nion, ns, nk, nlm, nb_win), dtype=np.complex128)
 
 # Select projectors for a given energy window
         ns_band = self.ib_win.shape[1]
@@ -288,11 +290,7 @@ class ProjectorShell:
                 ib2 = self.ib_win[ik, is_b, 1] + 1
                 ib1_win = ib1 - self.nb_min
                 ib2_win = ib2 - self.nb_min
-                self.proj_win[:, isp, ik, ib1_win:ib2_win, :] = self.proj_arr[:, isp, ik, ib1:ib2, :]
+                self.proj_win[:, isp, ik, :, ib1_win:ib2_win] = self.proj_arr[:, isp, ik, :, ib1:ib2]
-
-# !!! This sucks but I have to change the order of 'ib' and 'ilm' indices here
-# This should perhaps be done right after the projector array is read from PLOCAR
-        self.proj_win.transpose((0, 1, 2, 4, 3))
 
 
 def generate_ortho_plos(conf_pars, vasp_data):

python/converters/vasp/test/plotools/test_projgroups.py

@@ -17,39 +17,59 @@ class TestProjectorGroup(mytest.MyTestCase):
     ProjectorGroup(sh_pars, proj_raw)
 
     Scenarios:
-    - compare output for a correct input
+    - test output for a correct input
+    - test the output of 'orthogonalization()' (sanity check)
     """
+    def setUp(self):
+        conf_file = 'example.cfg'
+        self.pars = ConfigParameters(conf_file)
+        self.pars.parse_input()
+        self.vasp_data = vaspio.VaspData('./')
+
+        efermi = self.vasp_data.doscar.efermi
+        eigvals = self.vasp_data.eigenval.eigs - efermi
+
+        self.shells = [ProjectorShell(self.pars.shells[0], self.vasp_data.plocar.plo)]
+        self.proj_gr = ProjectorGroup(self.pars.groups[0], self.shells, eigvals)
+
 # Scenario 1
     def test_example(self):
-        conf_file = 'example.cfg'
-        pars = ConfigParameters(conf_file)
-        pars.parse_input()
-        print pars.groups
-        vasp_data = vaspio.VaspData('./')
-
-        efermi = vasp_data.doscar.efermi
-        eigvals = vasp_data.eigenval.eigs - efermi
-
-        shells = [ProjectorShell(pars.shells[0], vasp_data.plocar.plo)]
-        proj_gr = ProjectorGroup(pars.groups[0], shells, eigvals)
-
 #        proj_sh.select_projectors(ib_win, nb_min, nb_max)
 #
         testout = 'projgroups.out.test'
-        nion, ns, nk, nbtot, nlm = proj_gr.shells[0].proj_win.shape
+        nion, ns, nk, nlm, nbtot = self.proj_gr.shells[0].proj_win.shape
         with open(testout, 'wt') as f:
-            f.write("pars: %s\n"%(pars.groups[0]))
+            f.write("pars: %s\n"%(self.pars.groups[0]))
             for ion in xrange(nion):
                 for isp in xrange(ns):
                     for ik in xrange(nk):
-                        ib1 = proj_gr.ib_win[ik, 0, 0]
+                        ib1 = self.proj_gr.ib_win[ik, 0, 0]
-                        ib2 = proj_gr.ib_win[ik, 0, 1]
+                        ib2 = self.proj_gr.ib_win[ik, 0, 1]
                         f.write("%i  %i\n"%(ib1, ib2))
-                        for ib in xrange(ib2 - proj_gr.nb_min + 1):
+                        for ib in xrange(ib2 - self.proj_gr.nb_min + 1):
                             for ilm in xrange(nlm):
-                                p = proj_gr.shells[0].proj_win[ion, isp, ik, ib, ilm]
+                                p = self.proj_gr.shells[0].proj_win[ion, isp, ik, ilm, ib]
+                                f.write("%5i  %s\n"%(ilm+1, p))
+ 
+# Scenario 2
+    def test_ortho(self):
+        self.proj_gr.orthogonalize()
+
+        testout = 'projortho.out.test'
+        nion, ns, nk, nlm, nbtot = self.proj_gr.shells[0].proj_win.shape
+        with open(testout, 'wt') as f:
+            f.write("pars: %s\n"%(self.pars.groups[0]))
+            for ion in xrange(nion):
+                for isp in xrange(ns):
+                    for ik in xrange(nk):
+                        ib1 = self.proj_gr.ib_win[ik, 0, 0]
+                        ib2 = self.proj_gr.ib_win[ik, 0, 1]
+                        f.write("%i  %i\n"%(ib1, ib2))
+                        for ib in xrange(ib2 - self.proj_gr.nb_min + 1):
+                            for ilm in xrange(nlm):
+                                p = self.proj_gr.shells[0].proj_win[ion, isp, ik, ilm, ib]
                                 f.write("%5i  %s\n"%(ilm+1, p))
 
-        expected_file = 'projgroups.out'
+        expected_file = 'projortho.out'
         self.assertFileEqual(testout, expected_file)
- 
+

python/converters/vasp/test/plotools/test_projshells.py

@@ -37,7 +37,7 @@ class TestProjectorShell(mytest.MyTestCase):
         proj_sh.select_projectors(ib_win, nb_min, nb_max)
 
         testout = 'projshells.out.test'
-        nion, ns, nk, nbtot, nlm = proj_sh.proj_win.shape
+        nion, ns, nk, nlm, nbtot = proj_sh.proj_win.shape
         with open(testout, 'wt') as f:
             f.write("pars: %s\n"%(pars.shells[0]))
             for ion in xrange(nion):
@@ -48,7 +48,7 @@ class TestProjectorShell(mytest.MyTestCase):
                         f.write("%i  %i\n"%(ib1, ib2))
                         for ib in xrange(ib2 - nb_min + 1):
                             for ilm in xrange(nlm):
-                                p = proj_sh.proj_win[ion, isp, ik, ib, ilm]
+                                p = proj_sh.proj_win[ion, isp, ik, ilm, ib]
                                 f.write("%5i  %s\n"%(ilm+1, p))
 
         expected_file = 'projshells.out'