spectral routines plotpt3d update

python/triqs_dft_tools/converters/elktools/elk_converter_tools.py

@@ -297,8 +297,9 @@ class ElkConverterTools:
         return dy
         
     def plotpt3d(self,n_k,vkl,n_symm,symlat,grid3d,ngrid):
-        #import time
         import triqs.utility.mpi as mpi
+        #import time
+        #st = time.time()
         #default vector tolerance used in Elk. This should not be altered.
         epslat=1E-6
         tol=int(numpy.log10(1/epslat))
@@ -307,82 +308,54 @@ class ElkConverterTools:
         nk = ngrid[0]*ngrid[1]*ngrid[2]
         BZvkl = numpy.zeros([nk,3], float)
         BZvkl[:,:] = None
-        vklir = numpy.zeros([3], float)
-        vklir[:] = None
-        #array which maps the new vkl to the symmetrically equivalent interface self.vkl
+        #array which maps the new vkl to the symmetrically equivalent interface vkl
         iknr = numpy.zeros([nk], int)
-        v = numpy.zeros([3], float)
         nk_ = 0
-        ik = 0
         vklIBZ = [self.v3frac(vkl[ik,:],epslat) for ik in range(n_k)]
         vklIBZ = numpy.array(vklIBZ)
-        #print(vklIBZ)
-        #st = time.time()
-        #loop over the number of grid points for each reciprocal lattice
-        for i2 in range(ngrid[2]):
-          t2=float(i2)/float(ngrid[2])
-          for i1 in range(ngrid[1]):
-            t1=float(i1)/float(ngrid[1])
-            for i0 in range(ngrid[0]):
-              t0=float(i0)/float(ngrid[0])
-              #br = None
-              #calculate Brillouin zone lattice vector
-              v = t0*b[0,:]+t1*b[1,:]+t2*b[2,:]+grid3d[0,:]
-              BZvkl[ik,:] = v.copy()
-              ik += 1
-
-        #check if generated points are symmetrically equivalent to interfaced vkl
-        #Note that this loop is the bottle neck for this routine hence the 
-        #parallelisation.
-        ikarray = numpy.array(range(nk))
-        for ik in mpi.slice_array(ikarray):
-              #apply translation to reduce back to first Brillouin zone
-              v = self.v3frac(BZvkl[ik,:].copy(),epslat)
-              br = None
-
-              if v.round(tol).tolist() in vklIBZ.copy().round(tol).tolist():
-                #Find index of v in self.vkl
-                ikk = vkl.copy().round(tol).tolist().index(v.round(tol).tolist())
-                iknr[ik] = ikk
-                #ikir = numpy.append(ikir,ikk)
-                #check if v is a irreducible vector and tally these vectors
-                if v.round(tol).tolist() not in vklir.round(tol).tolist():
-                  nk_+=1
-                  vklir = numpy.vstack((vklir,v))
-                continue
-              #if v is not in interface set, see if it's symmetrically equivalent to
-              #a vector in self.vkl
-              for isym in range(n_symm):
-                if numpy.allclose(symlat[isym][:,:],numpy.eye(3)):
-                    continue
-                v_symm=numpy.matmul(symlat[isym][:,:].transpose(),v)
-                v_symm=self.v3frac(v_symm,epslat)
-                if v_symm.round(tol).tolist() in vklIBZ.copy().round(tol).tolist():
-                   ikk = vkl.copy().round(tol).tolist().index(v_symm.round(tol).tolist())
-                   #ikir = numpy.append(ikir,ikk)
-                   iknr[ik] = ikk
-                   br = 1
-                   break
-              if br == 1: continue
-              #if v is not symmetrically equivalent, then wrong input mesh.
-              mpi.report('No symmetrically equavilent vector in interface vkl set')
-              assert 0, "input grid does not generate interfaced reciprocal vectors"
-        #collect required variables and arrays (initialised to zero) from all threads.
-        nk_ = mpi.all_reduce(mpi.world, nk_, lambda x, y: x + y)
-        iknr = mpi.all_reduce(mpi.world, iknr, lambda x, y: x + y)
+               
+        #generate mesh grid      
+        i0, i1, i2 = numpy.meshgrid(numpy.arange(ngrid[0]), numpy.arange(ngrid[1]), 
+                                    numpy.arange(ngrid[2]), indexing='ij')
+        #convert to floats
+        t0 = i0.astype(float)/ngrid[0]
+        t1 = i1.astype(float)/ngrid[1]
+        t2 = i2.astype(float)/ngrid[2]
+        #Calculate Brillouin zone lattice vectors
+        BZvkl[:, 0] = (t0*b[0,0]+t1*b[1, 0]+t2*b[2, 0]+grid3d[0, 0]).flatten()
+        BZvkl[:, 1] = (t0*b[0,1]+t1*b[1, 1]+t2*b[2, 1]+grid3d[0, 1]).flatten()
+        BZvkl[:, 2] = (t0*b[0,2]+t1*b[1, 2]+t2*b[2, 2]+grid3d[0, 2]).flatten() 
+        #check k-point has equivalent point dft-interfaced k-point list (this is a bottle neck for performance)
+        for ik in range(nk):
+          br = None                
+          v1 = self.v3frac(BZvkl[ik,:], epslat)
+          #see if v1 is symmetrically equivalent to a vector in IBZvkl
+          for isym in range(n_symm):
+            v_symm=numpy.matmul(symlat[isym][:,:].transpose(),v1)
+            v_symm=self.v3frac(v_symm,epslat)
+            if v_symm.round(tol).tolist() in vklIBZ.round(tol).tolist():
+              iknr[ik] = vkl.round(tol).tolist().index(v_symm.round(tol).tolist())
+              #if identity symmetry operation was used, this v1 must be in the IBZ vector set
+              if numpy.allclose(symlat[isym][:,:],numpy.eye(3)):
+                nk_+=1
+              br = 1
+              break                
+          if br == 1: continue
+          #if v1 is not symmetrically equivalent, then wrong input mesh.
+          mpi.report('No identity symmetry operator or symmetrically equivalent vector in interface vkl set')
+          assert 0, "input grid does not generate interfaced reciprocal vectors"
 
         #check that all the vectors from the interface are in this list of vectors
         if(nk_!=n_k):
-          mpi.report('Incorrect number of irreducible vectors with respect to self.vkl ')
+          mpi.report('Incorrect number of irreducible vectors with respect to vkl ')
           mpi.report('%s!=%s'%(nk_,n_k))
           assert 0, "input grid does not generate interfaced reciprocal vectors"
         #et = time.time()
-        #mpi.report(et-st,nk,nk_)
-        #assert 0, ""
+        #mpi.report(et-st,nk)
         return BZvkl, iknr, nk
     
     def bzfoldout(self,n_k,vkl,n_symm,symlat):
-        import triqs.utility.mpi as mpi
+        #import triqs.utility.mpi as mpi
         epslat=1E-6
         tol=int(numpy.log10(1/epslat))
         #new temporary arrays for expanding irreducible Brillouin zone
@@ -393,20 +366,13 @@ class ElkConverterTools:
         vkl2[0,:,:] = vkl[:,:].copy()
         iknr2[0,:] = iknr[:].copy()
         #expand irreducible Brillouin zone
-        ikarray = numpy.array(range(n_k))
-        for ik in mpi.slice_array(ikarray):
+        for ik in range(n_k):
           for isym in range(n_symm):
             #find point in BZ by symmetry operation
             v=numpy.matmul(symlat[isym][:,:].transpose(),vkl[ik,:])
-            #shift back in to range [0,1) - Elk specific
-            #v[:]=self.v3frac(v,epslat)
             #alter temporary arrays
             vkl2[isym,ik,:] = v[:]
             iknr2[isym,ik] = ik
-        # Collect data from mpi (adding to elements with zeros):
-        vkl2 = mpi.all_reduce(mpi.world, vkl2, lambda x, y: x + y)
-        iknr2 = mpi.all_reduce(mpi.world, iknr2, lambda x, y: x + y)
-        mpi.barrier()
         #flatten arrays
         BZvkl = vkl2.reshape(n_k*n_symm,3)
         iknr = iknr2.reshape(n_k*n_symm)

test/python/elk/elk_spectralcontours_convert.py

@@ -13,6 +13,7 @@ testdir = cwd+'/elk_spectralcontours_convert'
 #change to test directory
 os.chdir(testdir)
 
+#default k-mesh
 Converter = ElkConverter(filename='SrVO3', repacking=True)
 Converter.hdf_file = 'elk_spectralcontours_convert.out.h5'
 Converter.convert_dft_input()
@@ -27,16 +28,26 @@ fs_elk = SK.spectral_contours(broadening=0.01, mesh=mesh, with_Sigma=False, with
 omega_elk = SK.spectral_contours(broadening=0.01, mesh=mesh, with_Sigma=False, with_dc=False, FS=False, proj_type='wann', save_to_file=False)
 omega_range_elk = SK.spectral_contours(broadening=0.01, mesh=mesh, plot_range=(-0.5,2), with_Sigma=False, with_dc=False, FS=False, proj_type='wann', save_to_file=False)
 
+#user specified k-mesh - has to be same as used in elk.in
+Converter = ElkConverter(filename='SrVO3', repacking=True)
+Converter.hdf_file = 'elk_spectralcontours_convert.out.h5'
+ngrid=np.array([10,10,1],np.int_)
+kgrid=np.array([[0.0,0.0,0.0],[1.0,0.0,0.0],[0.0,1.0,0.0],[0.0,0.0,1.0]],np.float_)
+Converter.convert_contours_input(kgrid=kgrid,ngrid=ngrid)
+SK2 = SumkDFTTools(hdf_file='elk_spectralcontours_convert.out.h5', use_dft_blocks=True)
+fs_elk_user = SK2.spectral_contours(broadening=0.01, mesh=mesh, with_Sigma=False, with_dc=False, FS=True, proj_type='wann', save_to_file=False)
 
 if mpi.is_master_node():
 
-    #with HDFArchive('elk_fermisurface_convert.ref.h5', 'a') as ar:
+    #with HDFArchive('elk_spectralcontours_convert.ref.h5', 'a') as ar:
     #    ar['fs_elk'] = fs_elk
+    #    ar['fs_elk_user'] = fs_elk_user
     #    ar['omega_elk'] = omega_elk
     #    ar['omega_range_elk'] = omega_range_elk
     #    ar['mesh'] = [omin,omax,oN]
     with HDFArchive('elk_spectralcontours_convert.out.h5', 'a') as ar:
         ar['fs_elk'] = fs_elk
+        ar['fs_elk_user'] = fs_elk_user
         ar['omega_elk'] = omega_elk
         ar['omega_range_elk'] = omega_range_elk
         ar['mesh'] = [omin,omax,oN]