Added calculation of density matrices for a shell

A method 'density_matrix()' for evaluating a density matrix of a given shell
has been added to class ProjectorShell. It requires an ElectronicStructure
object as an input an by default produces a site- and spin-diagonal
part of the density matrix using the Fermi-weights obtained directly from VASP.
Ideally, this density matrix should coincide with the one calculated
within VASP itself (inside the LDA+U module).

Corresponding sanity test has been added, which shows only that the
calculation does not crash. Real numerical tests are needed.

python/converters/vasp/python/elstruct.py

@@ -15,6 +15,7 @@ class ElectronicStructure:
     - *efermi* (float) : Fermi level read from DOSCAR
     - *proj_raw* (array[complex]) : raw projectors from PLOCAR
     - *eigvals* (array[float]) : KS eigenvalues
+    - *ferw* (array[float]) : Fermi weights from VASP
     - *kmesh* (dict) : parameters of the `k`-mesh
     - *structure* (dict) : parameters of the crystal structure
     - *symmetry* (dict) : paramters of symmetry
@@ -29,12 +30,12 @@ class ElectronicStructure:
 
         self.kmesh = {'nktot': self.nktot}
         self.kmesh['kpoints'] = vasp_data.kpoints.kpts
-        self.kmesh['kwights'] = vasp_data.eigenval.kwghts
+        self.kmesh['kweights'] = vasp_data.eigenval.kwghts
         try:
             self.kmesh['ntet'] = vasp_data.kpoints.ntet
             self.kmesh['itet'] = vasp_data.kpoints.itet
             self.kmesh['volt'] = vasp_data.kpoints.volt
-        except NameError:
+        except AttributeError:
             pass
 
 # Note that one should not subtract this Fermi level from eigenvalues
@@ -55,6 +56,11 @@ class ElectronicStructure:
 # [see ProjectorGroup.orthogonalization()]
         self.proj_raw = vasp_data.plocar.plo.transpose((0, 1, 2, 4, 3))
 
+# In fact, Fermi weights do not depend on ions
+# FIXME: restructure the data in PLOCAR to remove the redundant dependency
+#        of 'ferw' on ions
+        self.ferw = vasp_data.plocar.ferw[0, :, :, :]
+
 # Check that the number of atoms is the same in PLOCAR and POSCAR
         natom_plo = vasp_data.plocar.params['nion']
         assert natom_plo == self.natom, "PLOCAR is inconsistent with POSCAR (number of atoms)"

python/converters/vasp/python/plotools.py

@@ -1,4 +1,5 @@
 
+import itertools as it
 import numpy as np
 
 class Projector:
@@ -292,6 +293,39 @@ class ProjectorShell:
                 ib2_win = ib2 - self.nb_min
                 self.proj_win[:, isp, ik, :, ib1_win:ib2_win] = self.proj_arr[:, isp, ik, :, ib1:ib2]
 
+################################################################################
+#
+# select_projectors
+#
+################################################################################
+    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"
+
+        occ_mats = np.zeros((ns, nion, nlm, nlm), dtype=np.float64)
+
+        kweights = el_struct.kmesh['kweights']
+        occnums = el_struct.ferw
+        ib1 = self.nb_min
+        ib2 = self.nb_max + 1
+        for isp in xrange(ns):
+            for ik, weight, occ in it.izip(it.count(), kweights, occnums[isp, :, :]):
+                for io in xrange(nion):
+                    proj_k = self.proj_win[isp, io, ik, ...]
+                    occ_mats[isp, io, :, :] += np.dot(proj_k * occ[ib1:ib2], 
+                                                 proj_k.conj().T).real * weight
+
+#        if not symops is None:
+#            occ_mats = symmetrize_matrix_set(occ_mats, symops, ions, perm_map)
+
+        return occ_mats
+
+
 
 def generate_ortho_plos(conf_pars, vasp_data):
     """

python/converters/vasp/python/vaspio.py

@@ -54,7 +54,9 @@ class Plocar:
     Properties
     ----------
 
-    plo (numpy.array((nion, ns, nk, nb, nlmmax))) : raw projectors
+    - *plo* (numpy.array((nion, ns, nk, nb, nlmmax))) : raw projectors
+    - *params* (dict) : parameters read from PLOCAR
+    - *ferw* (array(nion, ns, nk, nb)) : Fermi weights from VASP
     """
 
     def from_file(self, vasp_dir='./', plocar_filename='PLOCAR'):
@@ -206,11 +208,11 @@ class Kpoints:
     Properties
     ----------
 
-        nktot (int) : total number of k-points in the IBZ
+    - nktot (int) : total number of k-points in the IBZ
-        kpts (numpy.array((nktot, 3), dtype=float)) : k-point vectors (fractional coordinates)
+    - kpts (numpy.array((nktot, 3), dtype=float)) : k-point vectors (fractional coordinates)
-        ntet (int) : total number of k-point tetrahedra
+    - ntet (int) : total number of k-point tetrahedra
-        itet (numpy.array((ntet, 5), dtype=float) : array of tetrahedra
+    - itet (numpy.array((ntet, 5), dtype=float) : array of tetrahedra
-        volt (float) : volume of a tetrahedron (the k-grid is assumed to
+    - volt (float) : volume of a tetrahedron (the k-grid is assumed to
           be uniform)
     """
 #

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

@@ -46,7 +46,9 @@ class TestProjectorGroup(mytest.MyTestCase):
                         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):
+                        ib1w = ib1 - self.proj_gr.nb_min
+                        ib2w = ib2 - self.proj_gr.nb_min + 1
+                        for ib in xrange(ib1w, ib2w):
                             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))
@@ -65,7 +67,9 @@ class TestProjectorGroup(mytest.MyTestCase):
                         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):
+                        ib1w = ib1 - self.proj_gr.nb_min
+                        ib2w = ib2 - self.proj_gr.nb_min + 1
+                        for ib in xrange(ib1w, ib2w):
                             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))

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

@@ -1,6 +1,7 @@
 
 import numpy as np
 import vaspio
+import elstruct
 from inpconf import ConfigParameters
 from plotools import select_bands, ProjectorShell
 import mytest
@@ -18,6 +19,7 @@ class TestProjectorShell(mytest.MyTestCase):
 
     Scenarios:
     - compare output for a correct input
+    - test density matrix
     """
 # Scenario 1
     def test_example(self):
@@ -54,3 +56,31 @@ class TestProjectorShell(mytest.MyTestCase):
         expected_file = 'projshells.out'
         self.assertFileEqual(testout, expected_file)
 
+# Scenario 2
+    def test_dens_mat(self):
+        conf_file = 'example.cfg'
+        pars = ConfigParameters(conf_file)
+        pars.parse_input()
+        vasp_data = vaspio.VaspData('./')
+        el_struct = elstruct.ElectronicStructure(vasp_data)
+
+        efermi = el_struct.efermi
+        eigvals = el_struct.eigvals - efermi
+        emin = pars.groups[0]['emin']
+        emax = pars.groups[0]['emax']
+        ib_win, nb_min, nb_max = select_bands(eigvals, emin, emax)
+
+        proj_sh = ProjectorShell(pars.shells[0], vasp_data.plocar.plo)
+
+        proj_sh.select_projectors(ib_win, nb_min, nb_max)
+
+        dens_mat = proj_sh.density_matrix(el_struct)
+        print dens_mat
+
+        testout = 'densmat.out.test'
+        with open(testout, 'wt') as f:
+            f.write("density matrix: %s\n"%(dens_mat))
+
+        expected_file = 'densmat.out'
+        self.assertFileEqual(testout, expected_file)
+