Added site-dependent phases to projectors

Now the projectors are defined in agreement with the formulas
in Amadon et al. Specifically, the phase exp(-i k Q) for site Q
is included.

python/converters/plovasp/elstruct.py

@@ -116,6 +116,20 @@ class ElectronicStructure:
 # 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)"
+        self.structure = {'a_brav': vasp_data.poscar.a_brav}
+        self.structure['nqtot'] = vasp_data.poscar.nq
+        self.structure['ntypes'] = vasp_data.poscar.ntypes
+        self.structure['nq_types'] = vasp_data.poscar.nions
+# Concatenate coordinates grouped by type into one array
+        self.structure['qcoords'] = np.vstack(vasp_data.poscar.q_types)
+        self.structure['type_of_ion'] = vasp_data.poscar.type_of_ion
+# FIXME: This can be removed if ion coordinates are stored in a continuous array
+## Construct a map to access coordinates by index
+#        self.structure['ion_index'] = []
+#        for isort, nq in enumerate(self.structure['nq_types']):
+#            for iq in xrange(nq):
+#                self.structure['ion_index'].append((isort, iq))
+
 
     def debug_density_matrix(self):
         """

python/converters/plovasp/plotools.py

@@ -109,7 +109,7 @@ def generate_plo(conf_pars, el_struct):
     print "  Generating %i shell%s..."%(nshell, '' if nshell == 1 else 's')
     pshells = []
     for sh_par in conf_pars.shells:
-        pshell = ProjectorShell(sh_par, proj_raw, el_struct.proj_params, el_struct.nc_flag)
+        pshell = ProjectorShell(sh_par, proj_raw, el_struct.proj_params, el_struct.kmesh, el_struct.structure, el_struct.nc_flag)
         print
         print "    Shell         : %s"%(pshell.user_index)
         print "    Orbital l     : %i"%(pshell.lorb)

python/converters/plovasp/proj_shell.py

@@ -68,7 +68,7 @@ class ProjectorShell:
     - proj_raw (numpy.array) : array of raw projectors
 
     """
-    def __init__(self, sh_pars, proj_raw, proj_params, nc_flag):
+    def __init__(self, sh_pars, proj_raw, proj_params, kmesh, structure, nc_flag):
         self.lorb = sh_pars['lshell']
         self.ion_list = sh_pars['ion_list']
         self.user_index = sh_pars['user_index']
@@ -84,7 +84,7 @@ class ProjectorShell:
         self.ndim = self.extract_tmatrices(sh_pars)
         self.nion = len(self.ion_list)
 
-        self.extract_projectors(proj_raw, proj_params)
+        self.extract_projectors(proj_raw, proj_params, kmesh, structure)
 
 ################################################################################
 #
@@ -203,7 +203,7 @@ class ProjectorShell:
 # extract_projectors
 #
 ################################################################################
-    def extract_projectors(self, proj_raw, proj_params):
+    def extract_projectors(self, proj_raw, proj_params, kmesh, structure):
         """
         Extracts projectors for the given shell.
         
@@ -223,13 +223,16 @@ class ProjectorShell:
             ndim = self.tmatrices.shape[1]
             self.proj_arr = np.zeros((nion, ns, nk, ndim, nb), dtype=np.complex128)
             for ik in xrange(nk):
+                kp = kmesh['kpoints'][ik]
                 for io, ion in enumerate(self.ion_list):
                     proj_k = np.zeros((ns, nlm, nb), dtype=np.complex128)
+                    qcoord = structure['qcoords'][ion]
+                    kphase = np.exp(-2.0j * np.pi * np.dot(kp, qcoord))
                     for m in xrange(nlm):
 # Here we search for the index of the projector with the given isite/l/m indices
                         for ip, par in enumerate(proj_params):
                             if par['isite'] - 1 == ion and par['l'] == self.lorb and par['m'] == m:
-                                proj_k[:, m, :] = proj_raw[ip, :, ik, :]
+                                proj_k[:, m, :] = proj_raw[ip, :, ik, :] * kphase
                                 break
                     for isp in xrange(ns):
                         self.proj_arr[io, isp, ik, :, :] = np.dot(self.tmatrices[io, :, :], proj_k[isp, :, :])
@@ -238,11 +241,15 @@ class ProjectorShell:
 # No transformation: just copy the projectors as they are
             self.proj_arr = np.zeros((nion, ns, nk, nlm, nb), dtype=np.complex128)
             for io, ion in enumerate(self.ion_list):
+                qcoord = structure['qcoords'][ion]
                 for m in xrange(nlm):
 # Here we search for the index of the projector with the given isite/l/m indices
                     for ip, par in enumerate(proj_params):
                         if par['isite'] - 1 == ion and par['l'] == self.lorb and par['m'] == m:
-                            self.proj_arr[io, :, :, m, :] = proj_raw[ip, :, :, :]
+                            for ik in xrange(nk):
+                                kp = kmesh['kpoints'][ik]
+                                kphase = np.exp(-2.0j * np.pi * np.dot(kp, qcoord))
+                                self.proj_arr[io, :, :, m, :] = proj_raw[ip, :, :, :] * kphase
                             break