diff --git a/python/triqs_dft_tools/converters/elktools/elk_converter_tools.py b/python/triqs_dft_tools/converters/elktools/elk_converter_tools.py index c9b6e71e..91bf92af 100644 --- a/python/triqs_dft_tools/converters/elktools/elk_converter_tools.py +++ b/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) diff --git a/test/python/elk/elk_spectralcontours_convert.py b/test/python/elk/elk_spectralcontours_convert.py index 8e81596a..cdf69a0d 100644 --- a/test/python/elk/elk_spectralcontours_convert.py +++ b/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] diff --git a/test/python/elk/elk_spectralcontours_convert/elk_spectralcontours_convert.ref.h5 b/test/python/elk/elk_spectralcontours_convert/elk_spectralcontours_convert.ref.h5 index 9d34e53e..22d591da 100644 Binary files a/test/python/elk/elk_spectralcontours_convert/elk_spectralcontours_convert.ref.h5 and b/test/python/elk/elk_spectralcontours_convert/elk_spectralcontours_convert.ref.h5 differ diff --git a/test/python/elk/elk_spectralcontours_convert/elk_spectralcontours_convert.ref_orig.h5 b/test/python/elk/elk_spectralcontours_convert/elk_spectralcontours_convert.ref_orig.h5 new file mode 100644 index 00000000..9d34e53e Binary files /dev/null and b/test/python/elk/elk_spectralcontours_convert/elk_spectralcontours_convert.ref_orig.h5 differ