plovasp: added documentation and tests

python/converters/plovasp/inpconf.py

@@ -221,8 +221,8 @@ class ConfigParameters:
         than the second one.
         """
         ftmp = map(int, par_str.split())
-        assert len(ftmp) == 2, "BWINDOW must be specified by exactly two ints"
+        assert len(ftmp) == 2, "BANDS must be specified by exactly two ints"
-        assert ftmp[0] < ftmp[1], "The first int in BWINDOW must be smaller than the second one"
+        assert ftmp[0] < ftmp[1], "The first int in BANDS must be smaller than the second one"
         return tuple(ftmp)
 
 ################################################################################

python/converters/plovasp/plotools.py

@@ -180,7 +180,7 @@ def generate_plo(conf_pars, el_struct):
 
             emesh = np.linspace(dos_emin, dos_emax, n_points)
             for ish in pgroup.ishells:
-                if  not isinstance(pshells[pgroup.ishells[ish]],ComplementShell):
+                if  not isinstance(pshells[pgroup.ishells[ish]],ComplementShell) or True:
                     print "  Shell %i"%(ish + 1)
                     dos = pshells[pgroup.ishells[ish]].density_of_states(el_struct, emesh)
                     de = emesh[1] - emesh[0]

python/converters/plovasp/proj_group.py

@@ -70,8 +70,8 @@ class ProjectorGroup:
         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[:,:,0] = gr_pars['bands'][0]-1
-            ib_win[:,:,1] = gr_pars['bands'][1]
+            ib_win[:,:,1] = gr_pars['bands'][1]-1
             ib_min = gr_pars['bands'][0]
             ib_max = gr_pars['bands'][1]
             
@@ -85,11 +85,10 @@ class ProjectorGroup:
         
         
         if gr_pars['complement']:
-            n_bands = self.ib_win[:,:,1] - self.ib_win[:,:,0]
+            n_bands = self.ib_win[:,:,1] - self.ib_win[:,:,0]+1
             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]:
-            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"
@@ -187,14 +186,17 @@ class ProjectorGroup:
 ################################################################################
     def calc_hk(self, eigvals):
         """
-        Calculate H(k) for a group by applying the projectors to the eigenvalues.
+        Calculate H(k) for a group by applying the projectors P
+        to the eigenvalues eps.
+        
+        H_ij(k) = sum_l P*_il eps_l P_lj
         
         """
 
         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
@@ -204,6 +206,7 @@ class ProjectorGroup:
                 bmax = self.ib_win[ik, isp, 1]+1
                 
                 nb = bmax - bmin 
+                p_mat = np.zeros((ndim, nb), dtype=np.complex128)
                 #print(bmin,bmax,nb)
 # Combine all projectors of the group to one block projector
                 for bl_map in block_maps:
@@ -214,14 +217,9 @@ class ProjectorGroup:
                         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],
+                self.hk[isp,ik,:,:] = np.dot(p_mat.conjugate()*eigvals[ik,bmin:bmax,isp],
-                                                p_mat.transpose().conjugate())
+                                        p_mat.transpose())
-                    #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())
 
 
 ################################################################################
@@ -232,6 +230,17 @@ class ProjectorGroup:
     def complement(self,eigvals):
         """
         Calculate the complement for a group of projectors.
+            
+        This leads to quadtratic projectors P = <l|n> by using a Gram-Schmidt.
+        
+        The projector on the orthogonal complement of the existing projectors
+        {|l>} is P^u = 1 - sum_l |l><l|
+        We get candidates for complement projectors by applying P^u to a Bloch
+        state |n>: |l*> = P^u |n>. For numerical stability we select that Bloch
+        state which leads to the |l*> with the largest norm (that corresponds to
+        that Bloch state with the smallest overlap with the space spanned by {|l>})
+        We normalize |l*> and add it to {|l>}. We do so untill we have as many
+        |l> states as we have {|n>} states.
         
         """
 
@@ -242,42 +251,47 @@ class ProjectorGroup:
         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 isp in xrange(ns):
-                for ik in xrange(nk):
+            for ik in xrange(nk):
-                    bmin = self.ib_win[ik, isp, 0]
+                bmin = self.ib_win[ik, isp, 0]
-                    bmax = self.ib_win[ik, isp, 1]+1
+                bmax = self.ib_win[ik, isp, 1]+1
-                    
+                
-                    nb = bmax - bmin 
+                nb = bmax - bmin 
-                    #print(bmin,bmax,nb)
+# Combine all projectors of the group to one block projector
-    # Combine all projectors of the group to one block projector
+                for bl_map in block_maps:
-                    for bl_map in block_maps:
+                    p_mat[:, :] = 0.0j  # !!! Clean-up from the last k-point and block!
-                        p_mat[:, :] = 0.0j  # !!! Clean-up from the last k-point and block!
+                    for ibl, block in enumerate(bl_map):
-                        for ibl, block in enumerate(bl_map):
+                        i1, i2 = block['bmat_range']
-                            i1, i2 = block['bmat_range']
+                        ish, ion = block['shell_ion']
-                            ish, ion = block['shell_ion']
+                        nlm = i2 - i1 + 1
-                            nlm = i2 - i1 + 1
+                        shell = self.shells[ish]
-                            shell = self.shells[ish]
+                        p_mat[i1:i2, :nb] = shell.proj_win[ion, isp, ik, :nlm, :nb]
-                            p_mat[i1:i2, :nb] = shell.proj_win[ion, isp, ik, :nlm, :nb]
+                orbs_done = 1*ndim
-                    p_full[0,isp,ik,:ndim,:] = p_mat
+                p_full[0,isp,ik,:ndim,:] = p_mat
+                while orbs_done < self.nb_max:                    
                     proj_work = p_full[0,isp,ik,:,:] 
+#We calculate the overlap of all bloch states: sum_l <n|l><l|m>
                     overlap = np.dot(proj_work.transpose().conjugate(),proj_work)
+# work is the projector onto the orthogonal complment <n| ( 1 - sum_l |l><l| ) |m>
                     work = np.eye(self.nb_max) - overlap
+# calculate the norm of the projected bloch function
                     norm = np.sqrt(np.sum(work*work.transpose().conjugate(),axis=1))
+# select the bloch function leading to the largest norm
                     max_ind = np.argmax(norm)
+# normalize and put it to the projectors
                     p_full[0,isp,ik,orbs_done,:] = work[max_ind]/norm[max_ind]
-                    
+                    orbs_done += 1
-            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
+        sh_pars['ib_min']  = bmin
+        sh_pars['ib_max']  = bmax
+        sh_pars['ib_win']  = self.ib_win
         self.shells.append(ComplementShell(sh_pars,p_full[:,:,:,ndim:,:],False))
         self.ishells.append(self.ishells[-1]+1)
                                                 

python/converters/plovasp/proj_shell.py

@@ -473,7 +473,11 @@ class ComplementShell(ProjectorShell):
         self.ions = sh_pars['ions']
         self.user_index = sh_pars['user_index']
         self.corr = sh_pars['corr']
-        self.nc_flag = nc_flag        
+        self.nc_flag = nc_flag   
+        
+        self.ib_min = sh_pars['ib_min']
+        self.ib_max = sh_pars['ib_max']
+        self.ib_win = sh_pars['ib_win']
 
 
         #self.lm1 = self.lorb**2
@@ -501,6 +505,6 @@ class ComplementShell(ProjectorShell):
     def density_matrix(self, el_struct, site_diag=True, spin_diag=True):
         raise Exception('not implemented')
 
-    def density_of_states(self, el_struct, emesh):
+    #def density_of_states(self, el_struct, emesh):
-        raise Exception('not implemented')
+    #    raise Exception('not implemented')
     

test/plovasp/inpconf/parse_shells_5.cfg

@@ -0,0 +1,13 @@
+[General]
+
+[Group 1]
+SHELLS = 1 2
+
+[Shell 1]
+LSHELL = 2
+IONS = 5..8
+
+[Shell 2]
+LSHELL = 1
+IONS = 1..4
+CORR = False

test/plovasp/inpconf/parse_shells_5.cfg~

@@ -0,0 +1,18 @@
+[General]
+
+[Group 1]
+SHELLS = 1 2
+
+[Shell 1]
+LSHELL = 2
+IONS = 5..8
+
+[Shell 2]
+LSHELL = 1
+IONS = 1..4
+
+TRANSFORM = 0.0  1.0  0.0
+            1.0  0.0  0.0
+            0.0  0.0  1.0
+
+CORR = False

test/plovasp/inpconf/test_general.py

@@ -30,7 +30,8 @@ class TestParseGeneral(arraytest.ArrayTestCase):
         conf_pars.parse_general()
         res = conf_pars.general
         expected = {'basename': 'test_base', 'efermi': 0.1,
-                    'dosmesh': {'n_points': 101, 'emin': -8.0, 'emax': 4.0}}
+                    'dosmesh': {'n_points': 101, 'emin': -8.0, 'emax': 4.0},
+                    'hk' : False}
         self.assertDictEqual(res, expected)
 
 

test/plovasp/inpconf/test_groups.py

@@ -39,9 +39,11 @@ class TestParseGroups(arraytest.ArrayTestCase):
         conf_pars.parse_groups()
         res = conf_pars.groups
         expected = [{'index': 1, 'shells': [1, 2], 'ewindow': (-7.6, 3.0),
-                     'normalize': True, 'normion': True},
+                     'normalize': True, 'normion': True,'complement': False},
                     {'index': 2, 'shells': [3], 'ewindow': (-1.6, 2.0),
-                     'normalize': True, 'normion': True}]
+                     'normalize': True, 'normion': True,'complement': False}]
+        print res
+        print expected
         self.assertListEqual(res, expected)
 
 

test/plovasp/inpconf/test_input.py

@@ -78,12 +78,12 @@ class TestParseInput(arraytest.ArrayTestCase):
         res = res.replace(" ","") # Remove spaces for comparison
 
         expected = r"""Shells:
-[{'ions':{'nion':4,'ion_list':[[4],[5],[6],[7]]},'user_index':1,'lshell':2},{'tmatrix':array([[0.,1.,0.],
+[{'ions':{'nion':4,'ion_list':[[4],[5],[6],[7]]},'user_index':1,'lshell':2,'corr':True},{'tmatrix':array([[0.,1.,0.],
 [1.,0.,0.],
-[0.,0.,1.]]),'ions':{'nion':4,'ion_list':[[0],[1],[2],[3]]},'user_index':2,'lshell':1},{'ions':{'nion':4,'ion_list':[[0],[1],[2],[3]]},'user_index':3,'lshell':3}]
+[0.,0.,1.]]),'ions':{'nion':4,'ion_list':[[0],[1],[2],[3]]},'user_index':2,'lshell':1,'corr':True},{'ions':{'nion':4,'ion_list':[[0],[1],[2],[3]]},'user_index':3,'lshell':3,'corr':True}]
 
 Groups:
-[{'normalize':True,'index':1,'ewindow':(-7.6,3.0),'normion':True,'shells':[0,1]},{'normalize':True,'index':2,'ewindow':(-1.6,2.0),'normion':True,'shells':[2]}]"""
+[{'normalize':True,'index':1,'ewindow':(-7.6,3.0),'shells':[0,1],'complement':False,'normion':True},{'normalize':True,'index':2,'ewindow':(-1.6,2.0),'shells':[2],'complement':False,'normion':True}]"""
 
         self.assertEqual(res, expected)
 
@@ -103,10 +103,10 @@ Groups:
         res = res.replace(" ","") # Remove spaces for comparison
 
         expected = r"""Shells:
-[{'ions':{'nion':4,'ion_list':[[4],[5],[6],[7]]},'user_index':1,'lshell':2}]
+[{'ions':{'nion':4,'ion_list':[[4],[5],[6],[7]]},'user_index':1,'lshell':2,'corr':True}]
 
 Groups:
-[{'normalize':True,'index':'1','ewindow':(-7.6,3.0),'shells':[0],'normion':True}]"""
+[{'normalize':True,'index':'1','ewindow':(-7.6,3.0),'normion':True,'complement':False,'shells':[0]}]"""
 
         self.assertEqual(res, expected)
 

test/plovasp/inpconf/test_shells.py

@@ -30,6 +30,8 @@ class TestParseShells(arraytest.ArrayTestCase):
       **raise** Exception
     - **if** two correct [Shell] sections are defined
       **return** a dictionary of shell parameters
+    - **if** two correct [Shell] sections (one has CORR=False are defined
+      **return** a dictionary of shell parameters
     """
 # Scenario 1
     def test_no_shell(self):
@@ -57,9 +59,9 @@ class TestParseShells(arraytest.ArrayTestCase):
         conf_pars = ConfigParameters(_rpath + 'parse_shells_4.cfg')
         conf_pars.parse_shells()
         res = conf_pars.shells
-        expected = [{'user_index': 1, 'lshell': 2, 'ions': {'nion': 4, 'ion_list': [[4],[5],[6],[7]]}},
+        expected = [{'user_index': 1, 'lshell': 2, 'ions': {'nion': 4, 'ion_list': [[4],[5],[6],[7]]},'corr': True},
                     {'user_index': 2, 'lshell': 1, 'ions': {'nion': 4, 'ion_list': [[0],[1],[2],[3]]},
-                        'tmatrix': np.array([[ 0.,  1.,  0.], [ 1.,  0.,  0.], [ 0.,  0.,  1.]])}]
+                        'tmatrix': np.array([[ 0.,  1.,  0.], [ 1.,  0.,  0.], [ 0.,  0.,  1.]]),'corr': True}]
 # ...lousy way to test equality of two dictionaries containing numpy arrays
         self.assertEqual(len(res), len(expected))
 
@@ -77,6 +79,24 @@ class TestParseShells(arraytest.ArrayTestCase):
 
         self.assertListEqual(res, expected)
 
+# Scenario 5
+    def test_two_shells_with_corr_false(self):
+        conf_pars = ConfigParameters(_rpath + 'parse_shells_5.cfg')
+        conf_pars.parse_shells()
+        res = conf_pars.shells
+        expected = [{'user_index': 1, 'lshell': 2, 'ions': {'nion': 4, 'ion_list': [[4],[5],[6],[7]]},'corr': True},
+                    {'user_index': 2, 'lshell': 1, 'ions': {'nion': 4, 'ion_list': [[0],[1],[2],[3]]},'corr': False}]
+        self.assertEqual(len(res), len(expected))
+
+        arr = res[0].pop('ions')
+        arr_exp = expected[0].pop('ions')
+        self.assertDictEqual(arr, arr_exp)
+
+        arr = res[1].pop('ions')
+        arr_exp = expected[1].pop('ions')
+        self.assertDictEqual(arr, arr_exp)
+
+        self.assertListEqual(res, expected)
 
 if __name__ == '__main__':
     import unittest

test/plovasp/inpconf/test_special_parsers.py

@@ -212,6 +212,59 @@ class TestParseEnergyWindow(arraytest.ArrayTestCase):
         with self.assertRaisesRegexp(AssertionError, err_mess):
             self.cpars.parse_energy_window('1.5 3.0 2.0')
 
+################################################################################
+#
+# TestParseBandWindow
+#
+################################################################################
+class TestParseBandWindow(arraytest.ArrayTestCase):
+    """
+    Function:
+
+    def parse_band_window(self, par_str)
+
+    Scenarios:
+
+    - **if** par_str == '1 10' **return** (1, 10)
+    - **if** par_str == '3.0 -1.5' **raise** AssertionError
+    - **if** par_str == '1.0' **raise** AssertionError
+    - **if** par_str == 'aaa' **raise** ValueError
+    - **if** par_str == '1.5 3.0 2.0' **raise** AssertionError
+    """
+    def setUp(self):
+        """
+        """
+# Dummy ConfigParameters object
+        self.cpars = ConfigParameters(_rpath + 'test1.cfg')
+
+# Scenario 1
+    def test_correct_range(self):
+        expected = (1, 10)
+        res = self.cpars.parse_band_window('1 10')
+        self.assertEqual(res, expected)
+
+# Scenario 2
+    def test_wrong_range(self):
+        err_mess = "The first int in BANDS"
+        with self.assertRaisesRegexp(AssertionError, err_mess):
+            self.cpars.parse_band_window('10 1')
+
+# Scenario 3
+    def test_one_float(self):
+        err_mess = "BANDS must be specified"
+        with self.assertRaisesRegexp(AssertionError, err_mess):
+            self.cpars.parse_band_window('1')
+
+# Scenario 4
+    def test_wrong_string(self):
+        with self.assertRaises(ValueError):
+            self.cpars.parse_band_window('aaa')
+
+# Scenario 5
+    def test_three_ints(self):
+        err_mess = "BANDS must be specified"
+        with self.assertRaisesRegexp(AssertionError, err_mess):
+            self.cpars.parse_band_window('1 2 3')
 
 ################################################################################
 #

test/plovasp/proj_shell/projshells.out

@@ -1,4 +1,4 @@
-pars: {'ions': {'nion': 1, 'ion_list': [[1]]}, 'user_index': 1, 'lshell': 2}
+pars: {'ions': {'nion': 1, 'ion_list': [[1]]}, 'user_index': 1, 'lshell': 2, 'corr': True}
 10  25
     1  0.000000  -0.000000
     2  0.000000  0.000000