Merge branch 'master' into aiida-parsing

Updated aiida-parsing branch with latest yambopy release

command_line/recipes.py

@@ -349,7 +349,7 @@ def add_qp(output,add=[],substract=[],addimg=[],verbose=False):
     addimgf=[f.name for f in addimg]
     filenames = addf+subf+addimgf
 
-    if len(filenames) is 0:
+    if len(filenames) == 0:
         raise ValueError('No files passed to function.')
 
 
@@ -686,4 +686,4 @@ def dimensions(array):
             else:
                 outVar[:] = varin[:]
 
-    fout.close()
+    fout.close()

qepy/pwxml.py

@@ -509,4 +509,121 @@ def spin_projection(self,spin_dir=3,folder='.',prefix='bands'):
                 ib1, ib2, ib3 = int(ib*10), int((ib+1)*10), int(ik*(nband/10+1)+2+ib)
                 self.spin_3[ik,ib1:ib2] = list( map(float,data_spin_3[ib3].split()))
 
+    def read_symmetries(self):
+        """
+        Read symmetry operations from data-file-schema.xml
+
+        Works for ibrav>0
+
+        NB: data-file-schema.xml has nrot and nsym with nsym<nrot.
+            They are stored in order with nsym first.
+
+        NB2: Most likely not working with symmorphic symmetries
+        """
+        #nsym
+        nrot = int(self.datafile_xml.findall("output/symmetries/nrot")[0].text.strip())        
+        self.nsym = int(self.datafile_xml.findall("output/symmetries/nsym")[0].text.strip())        
+        no_t_rev = ( self.datafile_xml.findall("input/symmetry_flags/no_t_rev")[0].text.strip() == "true" )
+        self.nsym_with_trev = 2*self.nsym 
+
+        # data
+        if not no_t_rev: sym_red = np.zeros((self.nsym_with_trev,3,3))
+        else:            sym_red = np.zeros((self.nsym,3,3))
+        symmetries = self.datafile_xml.findall("output/symmetries/symmetry/rotation") # All rotations
+        symmetries = symmetries[:self.nsym] # Retain only point group symmetries with no trev
+        #sym_info   = self.datafile_xml.findall("output/symmetries/symmetry/info")
+        #NB: sym_info[:].attrib['class'] contains irrep names if found
+
+        #read (non-symmorphic) symmetries
+        for i in range(self.nsym):
+            sym = np.array( [float(x) for x in symmetries[i].text.strip().split()] ).reshape(3,3)
+            sym_red[i] = sym.T # symmetries are saved as the transposed in the .xml 
+        self.sym_red = sym_red.astype(int)
+
+        #convert to c.c.
+        self.sym_car = self.sym_red_car()
+
+        #check for time reversal and apply it to sym_car
+        if not no_t_rev: self.apply_t_rev()
+
+    def apply_t_rev(self):
+        """
+        Add T*S rotation matrices
+        """
+        for n in range(self.nsym,self.nsym_with_trev): self.sym_car[n]=-1*self.sym_car[n-self.nsym]
 
+    def sym_red_car(self):
+        """
+        Transform symmetry ops. in Cartesian coordinates
+        """
+        lat = np.array(self.cell)
+        sym_car = np.zeros([len(self.sym_red),3,3],dtype=float)
+        for n,s in enumerate(self.sym_red):
+            sym_car[n] = np.dot( np.linalg.inv(lat), np.dot(s, lat ) ).T
+        return sym_car
+
+    def expand_kpoints_xml(self,atol=1e-6,expand_eigen=True,verbose=0):
+        """
+        Take a list of kpoints and symmetry operations and return the full brillouin zone
+        with the corresponding index in the irreducible brillouin zone
+
+        Expand also eigenvalues by default
+        """
+        alat      = np.array(self.acell)
+        kpts_ibz  = np.array(self.kpoints)
+        eigen_ibz = np.array(self.eigen1)
+        rlat      = np.array(self.rcell)
+
+        self.read_symmetries()
+
+        #check if the kpoints were already exapnded
+        kpoints_indices  = []
+        kpoints_full     = []
+        symmetry_indices = []
+
+        #kpoints in the full brillouin zone organized per index
+        kpoints_full_i = {}
+
+        #expand using symmetries
+        for nk,k in enumerate(kpts_ibz):
+            #if the index in not in the dictionary add a list
+            if nk not in kpoints_full_i:
+                kpoints_full_i[nk] = []
+
+            for ns,sym in enumerate(self.sym_car):
+
+                new_k = np.dot(sym,k)
+
+                #check if the point is inside the bounds
+                k_red = car_red([new_k],rlat)[0]
+                k_bz = (k_red+atol)%1
+
+                #if the vector is not in the list of this index add it
+                if not vec_in_list(k_bz,kpoints_full_i[nk]):
+                    kpoints_full_i[nk].append(k_bz)
+                    kpoints_full.append(new_k)
+                    kpoints_indices.append(nk)
+                    symmetry_indices.append(ns)
+                    continue
+
+        #calculate the weights of each of the kpoints in the irreducible brillouin zone
+        nkpoints_full = len(kpoints_full)
+        weights = np.zeros([nkpoints_full])
+        for nk in kpoints_full_i:
+            weights[nk] = float(len(kpoints_full_i[nk]))/nkpoints_full
+
+        if verbose: print("%d kpoints expanded to %d"%(self.nkpoints,len(kpoints_full)))
+
+        #set the variables
+        self.weights_ibz      = np.array(weights)
+        self.kpoints_indices  = np.array(kpoints_indices)
+        self.symmetry_indices = np.array(symmetry_indices)
+        self.nkbz             = nkpoints_full
+        #cartesian coordinates of QE
+        self.kpoints_bz       = np.array(kpoints_full)
+
+        if expand_eigen:
+
+            self.eigen_bz = np.zeros((self.nkbz,self.nbands))
+            for ik in range(self.nkbz): self.eigen_bz[ik,:] = eigen_ibz[self.kpoints_indices[ik],:]
+            if verbose: print("Eigenvalues also expanded.")

yambopy/__init__.py

@@ -20,6 +20,8 @@
         - YamboElectronsDB: read the electronic states from ns.db1
         - YamboQPDB: read the quasiparticle energies db ndb.QP
         - YamboGreenDB: read the green's functions calculated using yambo
+        - YamboExcitonsDB: read excitonic properties from a BSE calculation
+        - YamboKernelDB: read excitonic kernel in transition space
 
     bse
         - YamboExcitonWaveFunctionXSF: read the excitonic
@@ -57,6 +59,7 @@ class yambopyenv():
 from yambopy.dbs.excitondb import *
 from yambopy.dbs.wfdb import *
 from yambopy.dbs.elphondb import *
+from yambopy.dbs.bsekerneldb import *
 
 #input/output files
 from yambopy.io.inputfile import *

yambopy/dbs/electronsdb.py

@@ -262,20 +262,23 @@ def setFermiFixed(self,broad=1e-5):
 
     def energy_gaps(self,GWshift=0.):
         """
-        Calculate the enegy of the gap (by Fulvio Paleari)
+        Calculate the enegy of the gap and apply custom rigid shift
         """
         eiv = self.eigenvalues
         nv  = self.nbandsv
         nc  = self.nbandsc
 
+        # First apply shift if there is one
+        eiv[:,nv:]+=GWshift
+
+        # Then compute gaps
         homo = np.max(eiv[:,nv-1])
         lumo = np.min(eiv[:,nv])
         Egap = lumo-homo
         for k in eiv:
             if k[nv-1]==homo:
                 lumo_dir=k[nv]
         Edir = lumo_dir-homo
-        eiv[:,nv:]+=GWshift
 
         print('DFT Energy gap: %s eV'%Egap)
         print('DFT Direct gap: %s eV'%Edir)

yambopy/dbs/excitondb.py

@@ -1,16 +1,13 @@
-# Copyrigh (c) 2018, Henrique Miranda
-# All rights reserved.
-#
-# Thi file is part of the yambopy project
-#
-import os
-from itertools import product
-from yambopy import *
-from cmath import polar 
+
 from yambopy.units import *
-from yambopy.plot.plotting import add_fig_kwargs,BZ_hexagon
+from yambopy.plot.plotting import add_fig_kwargs,BZ_Wigner_Seitz
+from yambopy.plot.bandstructure import *
 from yambopy.lattice import replicate_red_kmesh, calculate_distances, get_path
 from yambopy.tools.funcs import gaussian, lorentzian
+from yambopy.dbs.savedb import *
+from yambopy.dbs.latticedb import *
+from yambopy.dbs.electronsdb import *
+
 
 class ExcitonList():
     """
@@ -54,7 +51,7 @@ class YamboExcitonDB(YamboSaveDB):
         Exciton eigenvectors are arranged as eigenvectors[i_exc, i_kvc]
         Transitions are unpacked in table[ i_k, i_v, i_c, i_s_c, i_s_v ] (last two are spin indices)
     """
-    def __init__(self,lattice,Qpt,eigenvalues,l_residual,r_residual,car_qpoint=None,q_cutoff=None,table=None,eigenvectors=None):
+    def __init__(self,lattice,Qpt,eigenvalues,l_residual,r_residual,spin_pol='no',car_qpoint=None,q_cutoff=None,table=None,eigenvectors=None):
         if not isinstance(lattice,YamboLatticeDB):
             raise ValueError('Invalid type for lattice argument. It must be YamboLatticeDB')
 
@@ -68,7 +65,7 @@ def __init__(self,lattice,Qpt,eigenvalues,l_residual,r_residual,car_qpoint=None,
         self.q_cutoff = q_cutoff
         self.table = table
         self.eigenvectors = eigenvectors
-        #self.spin_pol     = spin_pol
+        self.spin_pol = spin_pol
 
     @classmethod
     def from_db_file(cls,lattice,filename='ndb.BS_diago_Q1',folder='.'):
@@ -98,6 +95,7 @@ def from_db_file(cls,lattice,filename='ndb.BS_diago_Q1',folder='.'):
             if 'Q-point' in list(database.variables.keys()):
                 # Finite momentum
                 car_qpoint = database.variables['Q-point'][:]/lattice.alat
+            if Qpt=="1": car_qpoint = np.zeros(3)
 
             #energies
             eig =  database.variables['BS_Energies'][:]*ha2ev
@@ -114,7 +112,11 @@ def from_db_file(cls,lattice,filename='ndb.BS_diago_Q1',folder='.'):
 
             table = table
             eigenvectors = eigenvectors
-
+            spin_vars = [int(database.variables['SPIN_VARS'][:][0]), int(database.variables['SPIN_VARS'][:][1])]
+            if spin_vars[0] == 2 and spin_vars[1] == 1:
+               spin_pol = 'pol'
+            else:
+               spin_pol = 'no'
         # Check if Coulomb cutoff is present
         path_cutoff = os.path.join(path_filename.split('ndb',1)[0],'ndb.cutoff')  
         q_cutoff = None
@@ -124,7 +126,7 @@ def from_db_file(cls,lattice,filename='ndb.BS_diago_Q1',folder='.'):
                 bare_qpg = bare_qpg[:,:,0]+bare_qpg[:,:,1]*I
                 q_cutoff = np.abs(bare_qpg[0,int(Qpt)-1])
 
-        return cls(lattice,Qpt,eigenvalues,l_residual,r_residual,q_cutoff=q_cutoff,car_qpoint=car_qpoint,table=table,eigenvectors=eigenvectors)
+        return cls(lattice,Qpt,eigenvalues,l_residual,r_residual,spin_pol,q_cutoff=q_cutoff,car_qpoint=car_qpoint,table=table,eigenvectors=eigenvectors)
 
     @property
     def unique_vbands(self):
@@ -266,6 +268,7 @@ def exciton_bs(self,energies,path,excitons=(0,),debug=False):
             energies -> can be an instance of YamboSaveDB or YamboQBDB
             path     -> path in reduced coordinates in which to plot the band structure
             exciton  -> exciton index to plot
+            spin     -> ??
         """
         if self.eigenvectors is None:
             raise ValueError('This database does not contain Excitonic states,'
@@ -278,11 +281,12 @@ def exciton_bs(self,energies,path,excitons=(0,),debug=False):
 
         rep = list(range(-1,2))
         kpoints_rep, kpoints_idx_rep = replicate_red_kmesh(kpoints,repx=rep,repy=rep,repz=rep)
-        band_indexes = get_path(kpoints_rep,path)
+        band_indexes = get_path(kpoints_rep,path,debug=debug)
         band_kpoints = kpoints_rep[band_indexes] 
         band_indexes = kpoints_idx_rep[band_indexes]
 
         if debug:
+            import matplotlib.pyplot as plt
             for i,k in zip(band_indexes,band_kpoints):
                 x,y,z = k
                 plt.text(x,y,i) 
@@ -768,7 +772,7 @@ def plot_exciton_2D_ax(self,ax,excitons,f=None,mode='hexagon',limfactor=0.8,**kw
             excitons -> list of exciton indexes to plot
             f -> function to apply to the exciton weights. Ex. f=log will compute the 
                  log of th weight to enhance the small contributions
-            mode -> possible values are 'hexagon' to use hexagons as markers for the 
+            mode -> possible values are 'hexagon'/'square' to use hexagons/squares as markers for the 
                     weights plot and 'rbf' to interpolate the weights using radial basis functions.
             limfactor -> factor of the lattice parameter to choose the limits of the plot 
             scale -> size of the markers
@@ -780,13 +784,20 @@ def plot_exciton_2D_ax(self,ax,excitons,f=None,mode='hexagon',limfactor=0.8,**kw
         dlim = lim*1.1
         filtered_weights = [[xi,yi,di] for xi,yi,di in zip(x,y,weights_bz_sum) if -dlim<xi and xi<dlim and -dlim<yi and yi<dlim]
         x,y,weights_bz_sum = np.array(filtered_weights).T
+        # Add contours of BZ
+        ax.add_patch(BZ_Wigner_Seitz(self.lattice))
 
         #plotting
         if mode == 'hexagon': 
             scale = kwargs.pop('scale',1)
             ax.scatter(x,y,s=scale,marker='H',c=weights_bz_sum,rasterized=True,**kwargs)
             ax.set_xlim(-lim,lim)
             ax.set_ylim(-lim,lim)
+        elif mode == 'square': 
+            scale = kwargs.pop('scale',1)
+            ax.scatter(x,y,s=scale,marker='s',c=weights_bz_sum,rasterized=True,**kwargs)
+            ax.set_xlim(-lim,lim)
+            ax.set_ylim(-lim,lim)
         elif mode == 'rbf':
             from scipy.interpolate import Rbf
             npts = kwargs.pop('npts',100)
@@ -796,24 +807,16 @@ def plot_exciton_2D_ax(self,ax,excitons,f=None,mode='hexagon',limfactor=0.8,**kw
             weights_bz_sum = np.zeros([npts,npts])
             for col in range(npts):
                 weights_bz_sum[:,col] = rbfi(x,np.ones_like(x)*y[col])
-            ax.imshow(weights_bz_sum,interpolation=interp_method)
-
+            # NB we have to take the transpose of the imshow data to get the correct plot
+            ax.imshow(weights_bz_sum.T,interpolation=interp_method,extent=[-lim,lim,-lim,lim])
         title = kwargs.pop('title',str(excitons))
 
         ax.set_title(title)
         ax.set_aspect('equal')
         ax.set_xticks([])
         ax.set_yticks([])
-        ax.add_patch(BZ_hexagon(self.lattice.rlat))
+
         return ax
-
-    @add_fig_kwargs
-    def plot_exciton_2D(self,excitons,f=None,**kwargs):
-        import matplotlib.pyplot as plt
-        fig = plt.figure()
-        ax = fig.add_subplot(1,1,1)
-        self.plot_exciton_2D_ax(ax,excitons,f=f,**kwargs)
-        return fig
 
     def plot_nbrightest_2D(self,emin=0,emax=10,estep=0.001,broad=0.1,
                            mode='rbf',scale=3,nrows=2,ncols=2,eps=1e-5):
@@ -859,10 +862,13 @@ def get_exciton_bs(self,energies_db,path,excitons,size=1,space='bands',f=None,de
             raise ValueError('Path argument must be a instance of Path. Got %s instead'%type(path))
 
         if space == 'bands':
-            bands_kpoints, energies, weights = self.exciton_bs(energies_db, path.kpoints, excitons, debug)
-            nkpoints = len(bands_kpoints)
-            plot_energies = energies[:,self.start_band:self.mband]
-            plot_weights  = weights[:,self.start_band:self.mband]
+            if self.spin_pol=='no':
+               bands_kpoints, energies, weights = self.exciton_bs(energies_db, path.kpoints, excitons, debug)
+               nkpoints = len(bands_kpoints)
+               plot_energies = energies[:,self.start_band:self.mband]
+               plot_weights  = weights[:,self.start_band:self.mband]
+        #    elif spin_pol=='pol':
+
         else:
             raise NotImplementedError('TODO')
             eh_size = len(self.unique_vbands)*len(self.unique_cbands)
@@ -967,7 +973,6 @@ def interpolate(self,energies,path,excitons,lpratio=5,f=None,size=1,verbose=True
         for idx_bz,idx_ibz in enumerate(lattice.kpoints_indexes):
             ibz_weights[idx_ibz,:] = weights[idx_bz,:] 
             ibz_kpoints[idx_ibz] = lattice.red_kpoints[idx_bz]
-
         #get eigenvalues along the path
         if isinstance(energies,(YamboSaveDB,YamboElectronsDB)):
             #ibz_energies = energies.eigenvalues[:,self.start_band:self.mband] Old version
@@ -983,7 +988,6 @@ def interpolate(self,energies,path,excitons,lpratio=5,f=None,size=1,verbose=True
         skw = SkwInterpolator(lpratio,ibz_kpoints,ibz_energies[na,:,:],fermie,nelect,cell,symrel,time_rev,verbose=verbose)
         kpoints_path = path.get_klist()[:,:3]
         energies = skw.interp_kpts(kpoints_path).eigens
-
         #interpolate weights
         na = np.newaxis
         skw = SkwInterpolator(lpratio,ibz_kpoints,ibz_weights[na,:,:],fermie,nelect,cell,symrel,time_rev,verbose=verbose)
@@ -1059,7 +1063,6 @@ def interpolate_transitions(self,energies,path,excitons,lpratio=5,f=None,size=1,
 
         return exc_transitions
 
-
     def interpolate_spin(self,energies,spin_proj,path,excitons,lpratio=5,f=None,size=1,verbose=True,**kwargs):
         """ Interpolate exciton bandstructure using SKW interpolation from Abipy
         """

yambopy/dbs/qpdb.py

@@ -239,7 +239,10 @@ def interpolate(self,lattice,path,what='QP+KS',lpratio=5,valence=None,verbose=1,
         nelect = 0
         fermie = kwargs.pop('fermie',0)
 
-        #consistency check
+        #consistency check with lattice from YamboSaveDB
+        if len(lattice.kpts_iku)!=len(self.kpoints_iku):
+            print(len(lattice.kpts_iku),len(self.kpoints_iku))
+            raise ValueError("The QP database is not consistent with the lattice")
         if not np.isclose(lattice.kpts_iku,self.kpoints_iku).all():
             print(lattice.kpts_iku)
             print(self.kpoints_iku)