diff --git a/python/converters/hk_converter.py b/python/converters/hk_converter.py index 57f828de..f446632e 100644 --- a/python/converters/hk_converter.py +++ b/python/converters/hk_converter.py @@ -28,10 +28,10 @@ import string from math import sqrt -def Read_Fortran_File (filename): +def read_fortran_file (filename): """ Returns a generator that yields all numbers in the Fortran file as float, one by one""" import os.path - if not(os.path.exists(filename)) : raise IOError, "File %s does not exists"%filename + if not(os.path.exists(filename)) : raise IOError, "File %s does not exist."%filename for line in open(filename,'r') : for x in line.replace('D','E').replace('(',' ').replace(')',' ').replace(',',' ').split() : yield string.atof(x) @@ -40,22 +40,18 @@ def Read_Fortran_File (filename): class HkConverter: """ - Conversion from general H(k) file to an hdf5 file, that can be used as input for the SumK_LDA class. + Conversion from general H(k) file to an hdf5 file that can be used as input for the SumK_LDA class. """ def __init__(self, hk_file, hdf_file, lda_subgrp = 'SumK_LDA', symm_subgrp = 'SymmCorr', repacking = False): """ - Init of the class. Variable Filename gives the root of all filenames, e.g. case.ctqmcout, case.h5, and so + Init of the class. on. """ assert type(hk_file)==StringType,"hk_file must be a filename" self.hdf_file = hdf_file self.lda_file = hk_file - #self.Symm_file = Filename+'.symqmc' - #self.Parproj_file = Filename+'.parproj' - #self.Symmpar_file = Filename+'.sympar' - #self.Band_file = Filename+'.outband' self.lda_subgrp = lda_subgrp self.symm_subgrp = symm_subgrp @@ -72,12 +68,12 @@ class HkConverter: """ - if not (mpi.is_master_node()): return # do it only on master: + # Read and write only on the master node + if not (mpi.is_master_node()): return mpi.report("Reading input from %s..."%self.lda_file) - # Read and write only on Master!!! # R is a generator : each R.Next() will return the next number in the file - R = Read_Fortran_File(self.lda_file) + R = read_fortran_file(self.lda_file) try: energy_unit = 1.0 # the energy conversion factor is 1.0, we assume eV in files n_k = int(R.next()) # read the number of k points @@ -101,7 +97,6 @@ class HkConverter: self.inequiv_shells(corr_shells) # determine the number of inequivalent correlated shells, has to be known for further reading... - use_rotations = 0 rot_mat = [numpy.identity(corr_shells[icrsh][3],numpy.complex_) for icrsh in xrange(n_corr_shells)] rot_mat_time_inv = [0 for i in range(n_corr_shells)] @@ -109,16 +104,13 @@ class HkConverter: # Representative representations are read from file n_reps = [1 for i in range(self.n_inequiv_corr_shells)] dim_reps = [0 for i in range(self.n_inequiv_corr_shells)] - + T = [] for icrsh in range(self.n_inequiv_corr_shells): n_reps[icrsh] = int(R.next()) # number of representatives ("subsets"), e.g. t2g and eg dim_reps[icrsh] = [int(R.next()) for i in range(n_reps[icrsh])] # dimensions of the subsets - # The transformation matrix: - # it is of dimension 2l+1, it is taken to be standard d (as in Wien2k) - T = [] - for icrsh in range(self.n_inequiv_corr_shells): - #for ish in xrange(self.N_inequiv_corr_shells): + # The transformation matrix: + # is of dimension 2l+1, it is taken to be standard d (as in Wien2k) ll = 2*corr_shells[self.invshellmap[icrsh]][2]+1 lmax = ll * (corr_shells[self.invshellmap[icrsh]][4] + 1) T.append(numpy.zeros([lmax,lmax],numpy.complex_)) @@ -131,27 +123,21 @@ class HkConverter: # Spin blocks to be read: - n_spin_blocks = SP + 1 - SO # number of spins to read for Norbs and Ham, NOT Projectors + n_spin_blocs = SP + 1 - SO # number of spins to read for Norbs and Ham, NOT Projectors - # define the number of N_Orbitals for all k points: it is the number of total bands and independent of k! - n_orb = sum([ shells[ish][3] for ish in range(n_shells)]) - #n_orbitals = [ [n_orb for isp in range(n_spin_blocks)] for ik in xrange(n_k)] - n_orbitals = numpy.ones([n_k,n_spin_blocks],numpy.int) * n_orb - #print N_Orbitals + # define the number of n_orbitals for all k points: it is the number of total bands and independent of k! + n_orb = sum([ shells[ish][3] for ish in range(n_shells) ]) + n_orbitals = numpy.ones([n_k,n_spin_blocs],numpy.int) * n_orb # Initialise the projectors: - #proj_mat = [ [ [numpy.zeros([corr_shells[icrsh][3], n_orbitals[ik][isp]], numpy.complex_) - # for icrsh in range (n_corr_shells)] - # for isp in range(n_spin_blocks)] - # for ik in range(n_k) ] - proj_mat = numpy.zeros([n_k,n_spin_blocks,n_corr_shells,max(numpy.array(corr_shells)[:,3]),max(n_orbitals)],numpy.complex_) + proj_mat = numpy.zeros([n_k,n_spin_blocs,n_corr_shells,max(numpy.array(corr_shells)[:,3]),max(n_orbitals)],numpy.complex_) # Read the projectors from the file: for ik in xrange(n_k): for icrsh in range(n_corr_shells): - for isp in range(n_spin_blocks): + for isp in range(n_spin_blocs): # calculate the offset: offset = 0 @@ -169,9 +155,7 @@ class HkConverter: # now define the arrays for weights and hopping ... bz_weights = numpy.ones([n_k],numpy.float_)/ float(n_k) # w(k_index), default normalisation - #hopping = [ [numpy.zeros([n_orbitals[ik][isp],n_orbitals[ik][isp]],numpy.complex_) - # for isp in range(n_spin_blocks)] for ik in xrange(n_k) ] - hopping = numpy.zeros([n_k,n_spin_blocks,max(n_orbitals),max(n_orbitals)],numpy.complex_) + hopping = numpy.zeros([n_k,n_spin_blocs,max(n_orbitals),max(n_orbitals)],numpy.complex_) if (weights_in_file): # weights in the file @@ -181,11 +165,9 @@ class HkConverter: sm = sum(bz_weights) bz_weights[:] /= sm - # Grab the H - for ik in xrange(n_k) : - for isp in range(n_spin_blocks): - + for isp in range(n_spin_blocs): + for ik in xrange(n_k) : no = n_orbitals[ik,isp] if (first_real_part_matrix): @@ -220,24 +202,22 @@ class HkConverter: if ((only_upper_triangle)and(i!=j)): hopping[ik,isp,j,i] = hopping[ik,isp,i,j].conjugate() - #keep some things that we need for reading parproj: + # keep some things that we need for reading parproj: self.n_shells = n_shells self.shells = shells self.n_corr_shells = n_corr_shells self.corr_shells = corr_shells - self.n_spin_blocks = n_spin_blocks + self.n_spin_blocs = n_spin_blocs self.n_orbitals = n_orbitals self.n_k = n_k self.SO = SO self.SP = SP self.energy_unit = energy_unit except StopIteration : # a more explicit error if the file is corrupted. - raise "SumK_LDA : reading file HMLT_file failed!" + raise "HK Converter : reading file lda_file failed!" R.close() - #print Proj_Mat[0] - #----------------------------------------- # Store the input into HDF5: ar = HDFArchive(self.hdf_file,'a') @@ -276,7 +256,7 @@ class HkConverter: def __repack(self): """Calls the h5repack routine, in order to reduce the file size of the hdf5 archive. - Should only be used BEFORE the first invokation of HDF_Archive in the program, otherwise + Should only be used BEFORE the first invokation of HDFArchive in the program, otherwise the hdf5 linking is broken!!!""" import subprocess diff --git a/python/converters/wien2k_converter.py b/python/converters/wien2k_converter.py index 3c92b1f4..3d329112 100644 --- a/python/converters/wien2k_converter.py +++ b/python/converters/wien2k_converter.py @@ -30,7 +30,7 @@ import string def read_fortran_file (filename): """ Returns a generator that yields all numbers in the Fortran file as float, one by one""" import os.path - if not(os.path.exists(filename)) : raise IOError, "File %s does not exists"%filename + if not(os.path.exists(filename)) : raise IOError, "File %s does not exist."%filename for line in open(filename,'r') : for x in line.replace('D','E').split() : yield string.atof(x) @@ -39,13 +39,12 @@ def read_fortran_file (filename): class Wien2kConverter: """ - Conversion from Wien2k output to an hdf5 file, that can be used as input for the SumkLDA class. + Conversion from Wien2k output to an hdf5 file that can be used as input for the SumkLDA class. """ def __init__(self, filename, lda_subgrp = 'SumK_LDA', symm_subgrp = 'SymmCorr', repacking = False): """ - Init of the class. Variable filename gives the root of all filenames, e.g. case.ctqmcout, case.h5, and so - on. + Init of the class. Variable filename gives the root of all filenames, e.g. case.ctqmcout, case.h5, and so on. """ assert type(filename)==StringType,"LDA_file must be a filename" @@ -71,10 +70,10 @@ class Wien2kConverter: """ - if not (mpi.is_master_node()): return # do it only on master: + # Read and write only on the master node + if not (mpi.is_master_node()): return mpi.report("Reading input from %s..."%self.lda_file) - # Read and write only on Master!!! # R is a generator : each R.Next() will return the next number in the file R = read_fortran_file(self.lda_file) try: @@ -91,15 +90,13 @@ class Wien2kConverter: n_shells = int(R.next()) # number of shells (e.g. Fe d, As p, O p) in the unit cell, # corresponds to index R in formulas shells = [ [ int(R.next()) for i in range(4) ] for icrsh in range(n_shells) ] # reads iatom, sort, l, dim - #shells = numpy.array(shells) - + n_corr_shells = int(R.next()) # number of corr. shells (e.g. Fe d, Ce f) in the unit cell, # corresponds to index R in formulas # now read the information about the shells: corr_shells = [ [ int(R.next()) for i in range(6) ] for icrsh in range(n_corr_shells) ] # reads iatom, sort, l, dim, SO flag, irep self.inequiv_shells(corr_shells) # determine the number of inequivalent correlated shells, has to be known for further reading... - #corr_shells = numpy.array(corr_shells) use_rotations = 1 rot_mat = [numpy.identity(corr_shells[icrsh][3],numpy.complex_) for icrsh in xrange(n_corr_shells)] @@ -120,7 +117,7 @@ class Wien2kConverter: - # Read here the infos for the transformation of the basis: + # Read here the info for the transformation of the basis: n_reps = [1 for i in range(self.n_inequiv_corr_shells)] dim_reps = [0 for i in range(self.n_inequiv_corr_shells)] T = [] @@ -128,10 +125,8 @@ class Wien2kConverter: n_reps[icrsh] = int(R.next()) # number of representatives ("subsets"), e.g. t2g and eg dim_reps[icrsh] = [int(R.next()) for i in range(n_reps[icrsh])] # dimensions of the subsets - # The transformation matrix: - # it is of dimension 2l+1, if no SO, and 2*(2l+1) with SO!! - #T = [] - #for ish in xrange(self.n_inequiv_corr_shells): + # The transformation matrix: + # is of dimension 2l+1 without SO, and 2*(2l+1) with SO! ll = 2*corr_shells[self.invshellmap[icrsh]][2]+1 lmax = ll * (corr_shells[self.invshellmap[icrsh]][4] + 1) T.append(numpy.zeros([lmax,lmax],numpy.complex_)) @@ -151,21 +146,14 @@ class Wien2kConverter: # read the list of n_orbitals for all k points n_orbitals = numpy.zeros([n_k,n_spin_blocs],numpy.int) - #n_orbitals = [ [0 for isp in range(n_spin_blocs)] for ik in xrange(n_k)] for isp in range(n_spin_blocs): for ik in xrange(n_k): - #n_orbitals[ik][isp] = int(R.next()) n_orbitals[ik,isp] = int(R.next()) - #print n_orbitals # Initialise the projectors: - #proj_mat = [ [ [numpy.zeros([corr_shells[icrsh][3], n_orbitals[ik][isp]], numpy.complex_) - # for icrsh in range (n_corr_shells)] - # for isp in range(n_spin_blocs)] - # for ik in range(n_k) ] proj_mat = numpy.zeros([n_k,n_spin_blocs,n_corr_shells,max(numpy.array(corr_shells)[:,3]),max(n_orbitals)],numpy.complex_) - + # Read the projectors from the file: for ik in xrange(n_k): @@ -175,39 +163,34 @@ class Wien2kConverter: for isp in range(n_spin_blocs): for i in xrange(no): for j in xrange(n_orbitals[ik][isp]): - #proj_mat[ik][isp][icrsh][i,j] = R.next() proj_mat[ik,isp,icrsh,i,j] = R.next() # now Imag part: for isp in range(n_spin_blocs): for i in xrange(no): for j in xrange(n_orbitals[ik][isp]): - #proj_mat[ik][isp][icrsh][i,j] += 1j * R.next() proj_mat[ik,isp,icrsh,i,j] += 1j * R.next() - + # now define the arrays for weights and hopping ... bz_weights = numpy.ones([n_k],numpy.float_)/ float(n_k) # w(k_index), default normalisation - #hopping = [ [numpy.zeros([n_orbitals[ik][isp],n_orbitals[ik][isp]],numpy.complex_) - # for isp in range(n_spin_blocs)] for ik in xrange(n_k) ] hopping = numpy.zeros([n_k,n_spin_blocs,max(n_orbitals),max(n_orbitals)],numpy.complex_) - + # weights in the file for ik in xrange(n_k) : bz_weights[ik] = R.next() # if the sum over spins is in the weights, take it out again!! sm = sum(bz_weights) bz_weights[:] /= sm - + # Grab the H # we use now the convention of a DIAGONAL Hamiltonian!!!! for isp in range(n_spin_blocs): for ik in xrange(n_k) : - no = n_orbitals[ik][isp] + no = n_orbitals[ik,isp] for i in xrange(no): - #hopping[ik][isp][i,i] = R.next() * energy_unit hopping[ik,isp,i,i] = R.next() * energy_unit - #keep some things that we need for reading parproj: + # keep some things that we need for reading parproj: self.n_shells = n_shells self.shells = shells self.n_corr_shells = n_corr_shells @@ -219,12 +202,10 @@ class Wien2kConverter: self.SP = SP self.energy_unit = energy_unit except StopIteration : # a more explicit error if the file is corrupted. - raise "SumkLDA : reading file HMLT_file failed!" + raise "Wien2k_converter : reading file lda_file failed!" R.close() - #print proj_mat[0] - #----------------------------------------- # Store the input into HDF5: ar = HDFArchive(self.hdf_file,'a') @@ -279,25 +260,17 @@ class Wien2kConverter: for isp in range(self.n_spin_blocs) ] R = read_fortran_file(self.parproj_file) - #try: n_parproj = [int(R.next()) for i in range(self.n_shells)] n_parproj = numpy.array(n_parproj) # Initialise P, here a double list of matrices: - #proj_mat_pc = [ [ [ [numpy.zeros([self.shells[ish][3], self.n_orbitals[ik][isp]], numpy.complex_) - # for ir in range(n_parproj[ish])] - # for ish in range (self.n_shells) ] - # for isp in range(self.n_spin_blocs) ] - # for ik in range(self.n_k) ] - proj_mat_pc = numpy.zeros([self.n_k,self.n_spin_blocs,self.n_shells,max(n_parproj),max(numpy.array(self.shells)[:,3]),max(self.n_orbitals)],numpy.complex_) rot_mat_all = [numpy.identity(self.shells[ish][3],numpy.complex_) for ish in xrange(self.n_shells)] rot_mat_all_time_inv = [0 for i in range(self.n_shells)] for ish in range(self.n_shells): - #print ish # read first the projectors for this orbital: for ik in xrange(self.n_k): for ir in range(n_parproj[ish]): @@ -337,8 +310,6 @@ class Wien2kConverter: if (self.SP): rot_mat_all_time_inv[ish] = int(R.next()) - #except StopIteration : # a more explicit error if the file is corrupted. - # raise "Wien2kConverter: reading file for Projectors failed!" R.close() #----------------------------------------- @@ -378,10 +349,6 @@ class Wien2kConverter: n_orbitals[ik,isp] = int(R.next()) # Initialise the projectors: - #proj_mat = [ [ [numpy.zeros([self.corr_shells[icrsh][3], n_orbitals[ik][isp]], numpy.complex_) - # for icrsh in range (self.n_corr_shells)] - # for isp in range(self.n_spin_blocs)] - # for ik in range(n_k) ] proj_mat = numpy.zeros([n_k,self.n_spin_blocs,self.n_corr_shells,max(numpy.array(self.corr_shells)[:,3]),max(n_orbitals)],numpy.complex_) # Read the projectors from the file: @@ -399,8 +366,6 @@ class Wien2kConverter: for j in xrange(n_orbitals[ik,isp]): proj_mat[ik,isp,icrsh,i,j] += 1j * R.next() - #hopping = [ [numpy.zeros([n_orbitals[ik][isp],n_orbitals[ik][isp]],numpy.complex_) - # for isp in range(self.n_spin_blocs)] for ik in xrange(n_k) ] hopping = numpy.zeros([n_k,self.n_spin_blocs,max(n_orbitals),max(n_orbitals)],numpy.complex_) # Grab the H @@ -416,11 +381,6 @@ class Wien2kConverter: n_parproj = numpy.array(n_parproj) # Initialise P, here a double list of matrices: - #proj_mat_pc = [ [ [ [numpy.zeros([self.shells[ish][3], n_orbitals[ik][isp]], numpy.complex_) - # for ir in range(n_parproj[ish])] - # for ish in range (self.n_shells) ] - # for isp in range(self.n_spin_blocs) ] - # for ik in range(n_k) ] proj_mat_pc = numpy.zeros([n_k,self.n_spin_blocs,self.n_shells,max(n_parproj),max(numpy.array(self.shells)[:,3]),max(n_orbitals)],numpy.complex_) @@ -439,7 +399,7 @@ class Wien2kConverter: proj_mat_pc[ik,isp,ish,ir,i,j] += 1j * R.next() except StopIteration : # a more explicit error if the file is corrupted. - raise "SumkLDA : reading file HMLT_file failed!" + raise "Wien2k_converter : reading file band_file failed!" R.close() # reading done! @@ -448,18 +408,11 @@ class Wien2kConverter: # Store the input into HDF5: ar = HDFArchive(self.hdf_file,'a') if not (self.bands_subgrp in ar): ar.create_group(self.bands_subgrp) + # The subgroup containing the data. If it does not exist, it is created. # If it exists, the data is overwritten!!! thingstowrite = ['n_k','n_orbitals','proj_mat','hopping','n_parproj','proj_mat_pc'] for it in thingstowrite: exec "ar['%s']['%s'] = %s"%(self.bands_subgrp,it,it) - - #ar[self.bands_subgrp]['n_k'] = n_k - #ar[self.bands_subgrp]['n_orbitals'] = n_orbitals - #ar[self.bands_subgrp]['proj_mat'] = proj_mat - #self.proj_mat = proj_mat - #self.n_orbitals = n_orbitals - #self.n_k = n_k - #self.hopping = hopping del ar @@ -501,7 +454,7 @@ class Wien2kConverter: mat[in_s][orb][i,j] += 1j * R.next() # imaginary part # determine the inequivalent shells: - #SHOULD BE FINALLY REMOVED, PUT IT FOR ALL ORBITALS!!!!! + #SHOULD BE FINALLY REMOVED, PUT IT FOR ALL ORBITALS!!!!! (PS: FIXME?) #self.inequiv_shells(orbits) mat_tinv = [numpy.identity(orbits[orb][3],numpy.complex_) for orb in range(n_orbits)] @@ -519,7 +472,7 @@ class Wien2kConverter: except StopIteration : # a more explicit error if the file is corrupted. - raise "Symmetry : reading file failed!" + raise "Wien2k_converter : reading file symm_file failed!" R.close() diff --git a/python/sumk_lda.py b/python/sumk_lda.py index 4c11eed3..dd9950e3 100644 --- a/python/sumk_lda.py +++ b/python/sumk_lda.py @@ -21,18 +21,16 @@ ################################################################################ from types import * -#from pytriqs.applications.dft.symmetry import * from symmetry import * import numpy import pytriqs.utility.dichotomy as dichotomy from pytriqs.gf.local import * -#from pytriqs.applications.impurity_solvers.operators import * -from pytriqs.operators import * from pytriqs.archive import * import pytriqs.utility.mpi as mpi - from math import cos,sin +# FIXME PS: These two aren't used it seems +from pytriqs.operators.operators2 import * import string, pickle class SumkLDA: @@ -58,43 +56,40 @@ class SumkLDA: self.n_spin_blocks_gf = [2,1] self.Gupf = None self.h_field = h_field - + # read input from HDF: things_to_read = ['energy_unit','n_k','k_dep_projection','SP','SO','charge_below','density_required', 'symm_op','n_shells','shells','n_corr_shells','corr_shells','use_rotations','rot_mat', 'rot_mat_time_inv','n_reps','dim_reps','T','n_orbitals','proj_mat','bz_weights','hopping'] optional_things = ['gf_struct_solver','map_inv','map','chemical_potential','dc_imp','dc_energ','deg_shells'] - #ar=HDFArchive(self.hdf_file,'a') - #del ar - self.retval = self.read_input_from_hdf(subgrp=self.lda_data,things_to_read=things_to_read,optional_things=optional_things) - #ar=HDFArchive(self.hdf_file,'a') - #del ar - if (self.SO) and (abs(self.h_field)>0.000001): self.h_field=0.0 mpi.report("For SO, the external magnetic field is not implemented, setting it to 0!!") - - self.inequiv_shells(self.corr_shells) # determine the number of inequivalent correlated shells + + # determine the number of inequivalent correlated shells (self.n_inequiv_corr_shells) + # and related maps (self.shellmap, self.invshellmap) + self.inequiv_shells(self.corr_shells) # field to convert block_names to indices self.names_to_ind = [{}, {}] for ibl in range(2): - for inm in range(self.n_spin_blocks_gf[ibl]): + for inm in range(self.n_spin_blocks_gf[ibl]): self.names_to_ind[ibl][self.block_names[ibl][inm]] = inm * self.SP #(self.Nspinblocs-1) # GF structure used for the local things in the k sums - self.gf_struct_corr = [ [ (al, range( self.corr_shells[i][3])) for al in self.block_names[self.corr_shells[i][4]] ] + self.gf_struct_corr = [ [ (al, range( self.corr_shells[i][3])) for al in self.block_names[self.corr_shells[i][4]] ] for i in xrange(self.n_corr_shells) ] if not (self.retval['gf_struct_solver']): # No gf_struct was stored in HDF, so first set a standard one: - self.gf_struct_solver = [ [ (al, range( self.corr_shells[self.invshellmap[i]][3]) ) - for al in self.block_names[self.corr_shells[self.invshellmap[i]][4]] ] - for i in xrange(self.n_inequiv_corr_shells) ] + self.gf_struct_solver = [ dict([ (al, range(self.corr_shells[self.invshellmap[i]][3]) ) + for al in self.block_names[self.corr_shells[self.invshellmap[i]][4]] ]) + for i in range(self.n_inequiv_corr_shells) + ] self.map = [ {} for i in xrange(self.n_inequiv_corr_shells) ] self.map_inv = [ {} for i in xrange(self.n_inequiv_corr_shells) ] for i in xrange(self.n_inequiv_corr_shells): @@ -116,10 +111,10 @@ class SumkLDA: #mpi.report("Do the init for symm:") self.Symm_corr = Symmetry(hdf_file,subgroup=self.symm_corr_data) - # determine the smallest blocs, if wanted: + # Analyse the block structure and determine the smallest blocs, if desired if (use_lda_blocks): dm=self.analyse_BS() - + # now save things again to HDF5: if (mpi.is_master_node()): ar=HDFArchive(self.hdf_file,'a') @@ -127,21 +122,21 @@ class SumkLDA: del ar self.save() - - - + + + def read_input_from_hdf(self, subgrp, things_to_read, optional_things=[]): """ Reads data from the HDF file """ - + retval = True # init variables on all nodes: for it in things_to_read: exec "self.%s = 0"%it for it in optional_things: exec "self.%s = 0"%it - + if (mpi.is_master_node()): ar=HDFArchive(self.hdf_file,'a') if (subgrp in ar): @@ -152,7 +147,7 @@ class SumkLDA: else: mpi.report("Loading %s failed!"%it) retval = False - + if ((retval) and (len(optional_things)>0)): # if necessary things worked, now read optional things: retval = {} @@ -171,10 +166,10 @@ class SumkLDA: # now do the broadcasting: for it in things_to_read: exec "self.%s = mpi.bcast(self.%s)"%(it,it) for it in optional_things: exec "self.%s = mpi.bcast(self.%s)"%(it,it) - + retval = mpi.bcast(retval) - + return retval @@ -188,41 +183,41 @@ class SumkLDA: ar[self.lda_data]['chemical_potential'] = self.chemical_potential ar[self.lda_data]['dc_energ'] = self.dc_energ ar[self.lda_data]['dc_imp'] = self.dc_imp - del ar - + del ar + def load(self): """Loads some quantities from an HDF5 arxiv""" things_to_read=['chemical_potential','dc_imp','dc_energ'] - + retval = self.read_input_from_hdf(subgrp=self.lda_data,things_to_read=things_to_read) return retval def downfold(self,ik,icrsh,sig,gf_to_downfold,gf_inp): """Downfolding a block of the Greens function""" - + gf_downfolded = gf_inp.copy() isp = self.names_to_ind[self.SO][sig] # get spin index for proj. matrices dim = self.corr_shells[icrsh][3] n_orb = self.n_orbitals[ik,isp] - + gf_downfolded.from_L_G_R(self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb],gf_to_downfold,self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb].conjugate().transpose()) # downfolding G return gf_downfolded - + def upfold(self,ik,icrsh,sig,gf_to_upfold,gf_inp): """Upfolding a block of the Greens function""" gf_upfolded = gf_inp.copy() - + isp = self.names_to_ind[self.SO][sig] # get spin index for proj. matrices dim = self.corr_shells[icrsh][3] n_orb = self.n_orbitals[ik,isp] - - gf_upfolded.from_L_G_R(self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb].conjugate().transpose(),gf_to_upfold,self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb]) + + gf_upfolded.from_L_G_R(self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb].conjugate().transpose(),gf_to_upfold,self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb]) return gf_upfolded @@ -235,23 +230,23 @@ class SumkLDA: gf_rotated = gf_to_rotate.copy() if (direction=='toGlobal'): - #if (self.rot_mat_time_inv[icrsh]==1): gf_rotated <<= gf_rotated.transpose() - #gf_rotated.from_L_G_R(self.rot_mat[icrsh].transpose(),gf_rotated,self.rot_mat[icrsh].conjugate()) + if ((self.rot_mat_time_inv[icrsh]==1) and (self.SO)): gf_rotated <<= gf_rotated.transpose() gf_rotated.from_L_G_R(self.rot_mat[icrsh].conjugate(),gf_rotated,self.rot_mat[icrsh].transpose()) else: gf_rotated.from_L_G_R(self.rot_mat[icrsh],gf_rotated,self.rot_mat[icrsh].conjugate().transpose()) - + elif (direction=='toLocal'): + if ((self.rot_mat_time_inv[icrsh]==1)and(self.SO)): gf_rotated <<= gf_rotated.transpose() gf_rotated.from_L_G_R(self.rot_mat[icrsh].transpose(),gf_rotated,self.rot_mat[icrsh].conjugate()) else: gf_rotated.from_L_G_R(self.rot_mat[icrsh].conjugate().transpose(),gf_rotated,self.rot_mat[icrsh]) - + return gf_rotated - + def lattice_gf_matsubara(self,ik,mu,beta=40,with_Sigma=True): """Calculates the lattice Green function from the LDA hopping and the self energy at k-point number ik @@ -259,51 +254,77 @@ class SumkLDA: ntoi = self.names_to_ind[self.SO] bln = self.block_names[self.SO] - + if (not hasattr(self,"Sigma_imp")): with_Sigma=False - if (with_Sigma): + if (with_Sigma): stmp = self.add_dc() beta = self.Sigma_imp[0].mesh.beta #override beta if Sigma is present - - if (self.Gupf is None): - # first setting up of Gupf + +# +# if (self.Gupf is None): +# # first setting up of Gupf +# BS = [ range(self.n_orbitals[ik,ntoi[ib]]) for ib in bln ] +# gf_struct = [ (bln[ib], BS[ib]) for ib in range(self.n_spin_blocks_gf[self.SO]) ] +# a_list = [a for a,al in gf_struct] +# if (with_Sigma): #take the mesh from Sigma if necessary +# glist = lambda : [ GfImFreq(indices = al, mesh = self.Sigma_imp[0].mesh) for a,al in gf_struct] +# else: +# glist = lambda : [ GfImFreq(indices = al, beta = beta) for a,al in gf_struct] +# self.Gupf = BlockGf(name_list = a_list, block_list = glist(),make_copies=False) +# self.Gupf.zero() +# self.Gupf_id = self.Gupf.copy() +# self.Gupf_id <<= iOmega_n +# +# GFsize = [ gf.N1 for sig,gf in self.Gupf] +# unchangedsize = all( [ self.n_orbitals[ik,ntoi[bln[ib]]]==GFsize[ib] +# for ib in range(self.n_spin_blocks_gf[self.SO]) ] ) +# +# if ((not unchangedsize)or(self.Gupf.mesh.beta!=beta)): +# BS = [ range(self.n_orbitals[ik,ntoi[ib]]) for ib in bln ] +# gf_struct = [ (bln[ib], BS[ib]) for ib in range(self.n_spin_blocks_gf[self.SO]) ] +# a_list = [a for a,al in gf_struct] +# if (with_Sigma): +# glist = lambda : [ GfImFreq(indices = al, mesh = self.Sigma_imp[0].mesh) for a,al in gf_struct] +# else: +# glist = lambda : [ GfImFreq(indices = al, beta = beta) for a,al in gf_struct] +# self.Gupf = BlockGf(name_list = a_list, block_list = glist(),make_copies=False) +# self.Gupf.zero() +# self.Gupf_id = self.Gupf.copy() +# self.Gupf_id <<= iOmega_n +# +### +# FIXME PS Remove commented out code above if this works + # Do we need to set up Gupf? + set_up_Gupf = False + if self.Gupf == None: + set_up_Gupf = True + else: + GFsize = [ gf.N1 for sig,gf in self.Gupf] + unchangedsize = all( [ self.n_orbitals[ik,ntoi[bln[ib]]]==GFsize[ib] + for ib in range(self.n_spin_blocks_gf[self.SO]) ] ) + if ((not unchangedsize)or(self.Gupf.mesh.beta!=beta)): set_up_Gupf = True + + # Set up Gupf + if set_up_Gupf: BS = [ range(self.n_orbitals[ik,ntoi[ib]]) for ib in bln ] gf_struct = [ (bln[ib], BS[ib]) for ib in range(self.n_spin_blocks_gf[self.SO]) ] - a_list = [a for a,al in gf_struct] - if (with_Sigma): #take the mesh from Sigma if necessary - glist = lambda : [ GfImFreq(indices = al, mesh = self.Sigma_imp[0].mesh) for a,al in gf_struct] + a_list = [a for a,al in gf_struct] + if (with_Sigma): + glist = lambda : [ GfImFreq(indices = al, mesh = self.Sigma_imp[0].mesh) for a,al in gf_struct] else: glist = lambda : [ GfImFreq(indices = al, beta = beta) for a,al in gf_struct] self.Gupf = BlockGf(name_list = a_list, block_list = glist(),make_copies=False) self.Gupf.zero() self.Gupf_id = self.Gupf.copy() self.Gupf_id <<= iOmega_n +### - GFsize = [ gf.N1 for sig,gf in self.Gupf] - unchangedsize = all( [ self.n_orbitals[ik,ntoi[bln[ib]]]==GFsize[ib] - for ib in range(self.n_spin_blocks_gf[self.SO]) ] ) - - if ((not unchangedsize)or(self.Gupf.mesh.beta!=beta)): - BS = [ range(self.n_orbitals[ik,ntoi[ib]]) for ib in bln ] - gf_struct = [ (bln[ib], BS[ib]) for ib in range(self.n_spin_blocks_gf[self.SO]) ] - a_list = [a for a,al in gf_struct] - if (with_Sigma): - glist = lambda : [ GfImFreq(indices = al, mesh = self.Sigma_imp[0].mesh) for a,al in gf_struct] - else: - glist = lambda : [ GfImFreq(indices = al, beta = beta) for a,al in gf_struct] - self.Gupf = BlockGf(name_list = a_list, block_list = glist(),make_copies=False) - self.Gupf.zero() - self.Gupf_id = self.Gupf.copy() - self.Gupf_id <<= iOmega_n - - - idmat = [numpy.identity(self.n_orbitals[ik,ntoi[bl]],numpy.complex_) for bl in bln] - #for ibl in range(self.n_spin_blocks_gf[self.SO]): mupat[ibl] *= mu + idmat = [numpy.identity(self.n_orbitals[ik,ntoi[bl]],numpy.complex_) for bl in bln] self.Gupf <<= self.Gupf_id M = copy.deepcopy(idmat) - for ibl in range(self.n_spin_blocks_gf[self.SO]): + for ibl in range(self.n_spin_blocks_gf[self.SO]): ind = ntoi[bln[ibl]] n_orb = self.n_orbitals[ik,ind] M[ibl] = self.hopping[ik,ind,0:n_orb,0:n_orb] - (idmat[ibl]*mu) - (idmat[ibl] * self.h_field * (1-2*ibl)) @@ -320,11 +341,11 @@ class SumkLDA: def check_projectors(self): - dens_mat = [numpy.zeros([self.corr_shells[ish][3],self.corr_shells[ish][3]],numpy.complex_) + dens_mat = [numpy.zeros([self.corr_shells[ish][3],self.corr_shells[ish][3]],numpy.complex_) for ish in range(self.n_corr_shells)] - + for ik in range(self.n_k): - + for ish in range(self.n_corr_shells): dim = self.corr_shells[ish][3] n_orb = self.n_orbitals[ik,0] @@ -336,10 +357,10 @@ class SumkLDA: if (self.use_rotations): for icrsh in xrange(self.n_corr_shells): if (self.rot_mat_time_inv[icrsh]==1): dens_mat[icrsh] = dens_mat[icrsh].conjugate() - dens_mat[icrsh] = numpy.dot( numpy.dot(self.rot_mat[icrsh].conjugate().transpose(),dens_mat[icrsh]) , + dens_mat[icrsh] = numpy.dot( numpy.dot(self.rot_mat[icrsh].conjugate().transpose(),dens_mat[icrsh]) , self.rot_mat[icrsh] ) - - + + return dens_mat @@ -350,7 +371,7 @@ class SumkLDA: ntoi = self.names_to_ind[self.SO] bln = self.block_names[self.SO] - MMat = [numpy.zeros( [self.n_orbitals[0,ntoi[bl]],self.n_orbitals[0,ntoi[bl]]], numpy.complex_) for bl in bln] + MMat = [numpy.zeros( [self.n_orbitals[0,ntoi[bl]],self.n_orbitals[0,ntoi[bl]]], numpy.complex_) for bl in bln] dens_mat = [ {} for icrsh in xrange(self.n_corr_shells)] for icrsh in xrange(self.n_corr_shells): @@ -358,22 +379,22 @@ class SumkLDA: dens_mat[icrsh][bl] = numpy.zeros([self.corr_shells[icrsh][3],self.corr_shells[icrsh][3]], numpy.complex_) ikarray=numpy.array(range(self.n_k)) - + for ik in mpi.slice_array(ikarray): - - unchangedsize = all( [ self.n_orbitals[ik,ntoi[bln[ib]]]==len(MMat[ib]) + + unchangedsize = all( [ self.n_orbitals[ik,ntoi[bln[ib]]]==len(MMat[ib]) for ib in range(self.n_spin_blocks_gf[self.SO]) ] ) - + if (not unchangedsize): - MMat = [numpy.zeros( [self.n_orbitals[ik,ntoi[bl]],self.n_orbitals[ik,ntoi[bl]]], numpy.complex_) for bl in bln] + MMat = [numpy.zeros( [self.n_orbitals[ik,ntoi[bl]],self.n_orbitals[ik,ntoi[bl]]], numpy.complex_) for bl in bln] for ibl,bl in enumerate(bln): ind = ntoi[bl] for inu in range(self.n_orbitals[ik,ind]): - if ( (self.hopping[ik,ind,inu,inu]-self.h_field*(1-2*ibl)) < 0.0): + if ( (self.hopping[ik,ind,inu,inu]-self.h_field*(1-2*ibl)) < 0.0): MMat[ibl][inu,inu] = 1.0 else: - MMat[ibl][inu,inu] = 0.0 + MMat[ibl][inu,inu] = 0.0 for icrsh in range(self.n_corr_shells): @@ -381,9 +402,9 @@ class SumkLDA: isp = self.names_to_ind[self.corr_shells[icrsh][4]][bn] dim = self.corr_shells[icrsh][3] n_orb = self.n_orbitals[ik,isp] - + #print ik, bn, isp - dens_mat[icrsh][bn] += self.bz_weights[ik] * numpy.dot( numpy.dot(self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb],MMat[ibn]) , + dens_mat[icrsh][bn] += self.bz_weights[ik] * numpy.dot( numpy.dot(self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb],MMat[ibn]) , self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb].transpose().conjugate() ) # get data from nodes: @@ -392,7 +413,7 @@ class SumkLDA: dens_mat[icrsh][sig] = mpi.all_reduce(mpi.world,dens_mat[icrsh][sig],lambda x,y : x+y) mpi.barrier() - + if (self.symm_op!=0): dens_mat = self.Symm_corr.symmetrize(dens_mat) # Rotate to local coordinate system: @@ -400,15 +421,15 @@ class SumkLDA: for icrsh in xrange(self.n_corr_shells): for bn in dens_mat[icrsh]: if (self.rot_mat_time_inv[icrsh]==1): dens_mat[icrsh][bn] = dens_mat[icrsh][bn].conjugate() - dens_mat[icrsh][bn] = numpy.dot( numpy.dot(self.rot_mat[icrsh].conjugate().transpose(),dens_mat[icrsh][bn]) , + dens_mat[icrsh][bn] = numpy.dot( numpy.dot(self.rot_mat[icrsh].conjugate().transpose(),dens_mat[icrsh][bn]) , self.rot_mat[icrsh]) - + return dens_mat def density_gf(self,beta): - """Calculates the density without setting up Gloc. It is useful for Hubbard I, and very fast.""" + """Calculates the density without setting up Gloc. It is useful for Hubbard I, and very fast.""" dens_mat = [ {} for icrsh in xrange(self.n_corr_shells)] for icrsh in xrange(self.n_corr_shells): @@ -418,12 +439,12 @@ class SumkLDA: ikarray=numpy.array(range(self.n_k)) for ik in mpi.slice_array(ikarray): - + Gupf = self.lattice_gf_matsubara(ik=ik, beta=beta, mu=self.chemical_potential) Gupf *= self.bz_weights[ik] dm = Gupf.density() MMat = [dm[bl] for bl in self.block_names[self.SO]] - + for icrsh in range(self.n_corr_shells): for ibn,bn in enumerate(self.block_names[self.corr_shells[icrsh][4]]): isp = self.names_to_ind[self.corr_shells[icrsh][4]][bn] @@ -447,7 +468,7 @@ class SumkLDA: for icrsh in xrange(self.n_corr_shells): for bn in dens_mat[icrsh]: if (self.rot_mat_time_inv[icrsh]==1): dens_mat[icrsh][bn] = dens_mat[icrsh][bn].conjugate() - dens_mat[icrsh][bn] = numpy.dot( numpy.dot(self.rot_mat[icrsh].conjugate().transpose(),dens_mat[icrsh][bn]) , + dens_mat[icrsh][bn] = numpy.dot( numpy.dot(self.rot_mat[icrsh].conjugate().transpose(),dens_mat[icrsh][bn]) , self.rot_mat[icrsh] ) return dens_mat @@ -455,23 +476,22 @@ class SumkLDA: def analyse_BS(self, threshold = 0.00001, include_shells = None, dm = None): - """ Determines the Greens function block structure from the simple point integration""" + """ Determines the Green function block structure from simple point integration.""" if (dm==None): dm = self.simple_point_dens_mat() - + dens_mat = [dm[self.invshellmap[ish]] for ish in xrange(self.n_inequiv_corr_shells) ] - + if include_shells is None: include_shells=range(self.n_inequiv_corr_shells) for ish in include_shells: - #self.gf_struct_solver.append([]) - self.gf_struct_solver[ish] = [] + self.gf_struct_solver[ish] = {} gf_struct_temp = [] - + a_list = [a for a,al in self.gf_struct_corr[self.invshellmap[ish]] ] for a in a_list: - - dm = dens_mat[ish][a] + + dm = dens_mat[ish][a] dmbool = (abs(dm) > threshold) # gives an index list of entries larger that threshold offdiag = [] @@ -496,16 +516,16 @@ class SumkLDA: for i in range(NBlocs): blocs[i].sort() - self.gf_struct_solver[ish].append( ('%s%s'%(a,i),range(len(blocs[i]))) ) - gf_struct_temp.append( ('%s%s'%(a,i),blocs[i]) ) - - + self.gf_struct_solver[ish].update( [('%s_%s'%(a,i),range(len(blocs[i])))] ) + gf_struct_temp.append( ('%s_%s'%(a,i),blocs[i]) ) + + # map is the mapping of the blocs from the SK blocs to the CTQMC blocs: self.map[ish][a] = range(len(dmbool)) for ibl in range(NBlocs): for j in range(len(blocs[ibl])): - self.map[ish][a][blocs[ibl][j]] = '%s%s'%(a,ibl) - self.map_inv[ish]['%s%s'%(a,ibl)] = a + self.map[ish][a][blocs[ibl][j]] = '%s_%s'%(a,ibl) + self.map_inv[ish]['%s_%s'%(a,ibl)] = a # now calculate degeneracies of orbitals: @@ -546,13 +566,13 @@ class SumkLDA: except: mpi.report("deg_shells not stored, degeneracies not found") del ar - + return dens_mat - + def symm_deg_gf(self,gf_to_symm,orb): """Symmetrises a GF for the given degenerate shells self.deg_shells""" - + for degsh in self.deg_shells[orb]: #loop over degenerate shells: ss = gf_to_symm[degsh[0]].copy() @@ -560,7 +580,7 @@ class SumkLDA: Ndeg = len(degsh) for bl in degsh: ss += gf_to_symm[bl] / (1.0*Ndeg) for bl in degsh: gf_to_symm[bl] <<= ss - + def eff_atomic_levels(self): """Calculates the effective atomic levels needed as input for the Hubbard I Solver.""" @@ -570,14 +590,14 @@ class SumkLDA: for ish in range(self.n_inequiv_corr_shells): for bn in self.block_names[self.corr_shells[self.invshellmap[ish]][4]]: eff_atlevels[ish][bn] = numpy.identity(self.corr_shells[self.invshellmap[ish]][3], numpy.complex_) - + # Chemical Potential: - for ish in xrange(self.n_inequiv_corr_shells): + for ish in xrange(self.n_inequiv_corr_shells): for ii in eff_atlevels[ish]: eff_atlevels[ish][ii] *= -self.chemical_potential - + # double counting term: #if hasattr(self,"dc_imp"): - for ish in xrange(self.n_inequiv_corr_shells): + for ish in xrange(self.n_inequiv_corr_shells): for ii in eff_atlevels[ish]: eff_atlevels[ish][ii] -= self.dc_imp[self.invshellmap[ish]][ii] @@ -586,10 +606,10 @@ class SumkLDA: # calculate the sum over k. Does not depend on mu, so do it only once: self.Hsumk = [ {} for ish in range(self.n_corr_shells) ] for icrsh in range(self.n_corr_shells): - for bn in self.block_names[self.corr_shells[icrsh][4]]: + for bn in self.block_names[self.corr_shells[icrsh][4]]: dim = self.corr_shells[icrsh][3] #*(1+self.corr_shells[icrsh][4]) self.Hsumk[icrsh][bn] = numpy.zeros([dim,dim],numpy.complex_) - + for icrsh in range(self.n_corr_shells): dim = self.corr_shells[icrsh][3] for ibn, bn in enumerate(self.block_names[self.corr_shells[icrsh][4]]): @@ -598,7 +618,7 @@ class SumkLDA: n_orb = self.n_orbitals[ik,isp] MMat = numpy.identity(n_orb, numpy.complex_) MMat = self.hopping[ik,isp,0:n_orb,0:n_orb] - (1-2*ibn) * self.h_field * MMat - self.Hsumk[icrsh][bn] += self.bz_weights[ik] * numpy.dot( numpy.dot(self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb],MMat), #self.hopping[ik][isp]) , + self.Hsumk[icrsh][bn] += self.bz_weights[ik] * numpy.dot( numpy.dot(self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb],MMat), #self.hopping[ik][isp]) , self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb].conjugate().transpose() ) # symmetrisation: @@ -611,12 +631,12 @@ class SumkLDA: if (self.rot_mat_time_inv[icrsh]==1): self.Hsumk[icrsh][bn] = self.Hsumk[icrsh][bn].conjugate() #if (self.corr_shells[icrsh][4]==0): self.Hsumk[icrsh][bn] = self.Hsumk[icrsh][bn].conjugate() - - self.Hsumk[icrsh][bn] = numpy.dot( numpy.dot(self.rot_mat[icrsh].conjugate().transpose(),self.Hsumk[icrsh][bn]) , + + self.Hsumk[icrsh][bn] = numpy.dot( numpy.dot(self.rot_mat[icrsh].conjugate().transpose(),self.Hsumk[icrsh][bn]) , self.rot_mat[icrsh] ) - + # add to matrix: - for ish in xrange(self.n_inequiv_corr_shells): + for ish in xrange(self.n_inequiv_corr_shells): for bn in eff_atlevels[ish]: eff_atlevels[ish][bn] += self.Hsumk[self.invshellmap[ish]][bn] @@ -626,7 +646,7 @@ class SumkLDA: def __init_dc(self): - + # construct the density matrix dm_imp and double counting arrays #self.dm_imp = [ {} for i in xrange(self.n_corr_shells)] self.dc_imp = [ {} for i in xrange(self.n_corr_shells)] @@ -638,29 +658,12 @@ class SumkLDA: - def set_lichtenstein_dc(self,Sigma_imp): - """Sets a double counting term according to Lichtenstein et al. PRL2001""" - - assert 0,"Lichtenstein DC not supported any more!" - - - - def set_dc(self,dens_mat,U_interact,J_hund,orb=0,use_dc_formula=0,use_val=None): - """Sets the double counting term for inequiv orbital orb - use_dc_formula=0: LDA+U FLL double counting, use_dc_formula=1: Held's formula. - use_dc_formula=2: AMF - Be sure that you use the correct interaction Hamiltonian!""" - - - #if (not hasattr(self,"dc_imp")): self.__init_dc() - - - #dm = [ {} for i in xrange(self.n_corr_shells)] - #for i in xrange(self.n_corr_shells): - # l = self.corr_shells[i][3] #*(1+self.corr_shells[i][4]) - # for j in xrange(len(self.gf_struct_corr[i])): - # dm[i]['%s'%self.gf_struct_corr[i][j][0]] = numpy.zeros([l,l],numpy.float_) + """Sets the double counting term for inequiv orbital orb: + use_dc_formula=0: LDA+U FLL double counting, + use_dc_formula=1: Held's formula, + use_dc_formula=2: AMF. + Be sure that you are using the correct interaction Hamiltonian!""" for icrsh in xrange(self.n_corr_shells): @@ -671,21 +674,18 @@ class SumkLDA: # do this orbital Ncr = {} l = self.corr_shells[icrsh][3] #*(1+self.corr_shells[icrsh][4]) - + for j in xrange(len(self.gf_struct_corr[icrsh])): self.dc_imp[icrsh]['%s'%self.gf_struct_corr[icrsh][j][0]] = numpy.identity(l,numpy.float_) blname = self.gf_struct_corr[icrsh][j][0] Ncr[blname] = 0.0 - for a,al in self.gf_struct_solver[iorb]: - #for bl in self.map[iorb][blname]: + for a,al in self.gf_struct_solver[iorb].iteritems(): bl = self.map_inv[iorb][a] - #print 'bl, valiue = ',bl,dens_mat[a].real.trace() Ncr[bl] += dens_mat[a].real.trace() - #print 'Ncr=',Ncr M = self.corr_shells[icrsh][3] - + Ncrtot = 0.0 a_list = [a for a,al in self.gf_struct_corr[icrsh]] for bl in a_list: @@ -701,41 +701,37 @@ class SumkLDA: Ncr[bl] = Ncrtot / 2.0 if (use_val is None): - - if (use_dc_formula==0): + + if (use_dc_formula==0): # FLL self.dc_energ[icrsh] = U_interact / 2.0 * Ncrtot * (Ncrtot-1.0) for bl in a_list: Uav = U_interact*(Ncrtot-0.5) - J_hund*(Ncr[bl] - 0.5) - self.dc_imp[icrsh][bl] *= Uav + self.dc_imp[icrsh][bl] *= Uav self.dc_energ[icrsh] -= J_hund / 2.0 * (Ncr[bl]) * (Ncr[bl]-1.0) mpi.report("DC for shell %(icrsh)i and block %(bl)s = %(Uav)f"%locals()) - elif (use_dc_formula==1): + elif (use_dc_formula==1): # Held's formula, with U_interact the interorbital onsite interaction self.dc_energ[icrsh] = (U_interact + (M-1)*(U_interact-2.0*J_hund) + (M-1)*(U_interact-3.0*J_hund))/(2*M-1) / 2.0 * Ncrtot * (Ncrtot-1.0) - #self.dc_energ[icrsh] = (U_interact + J_hund * (2.0-(M-1)) / (2*M-1) ) / 2.0 * Ncrtot * (Ncrtot-1.0) for bl in a_list: - # Held's formula, with U_interact the interorbital onsite interaction Uav =(U_interact + (M-1)*(U_interact-2.0*J_hund) + (M-1)*(U_interact-3.0*J_hund))/(2*M-1) * (Ncrtot-0.5) - #Uav = (U_interact + J_hund * (2.0-(M-1)) / (2*M-1) ) * (Ncrtot-0.5) - self.dc_imp[icrsh][bl] *= Uav + self.dc_imp[icrsh][bl] *= Uav mpi.report("DC for shell %(icrsh)i and block %(bl)s = %(Uav)f"%locals()) - elif (use_dc_formula==2): + elif (use_dc_formula==2): # AMF self.dc_energ[icrsh] = 0.5 * U_interact * Ncrtot * Ncrtot for bl in a_list: - # AMF Uav = U_interact*(Ncrtot - Ncr[bl]/M) - J_hund * (Ncr[bl] - Ncr[bl]/M) self.dc_imp[icrsh][bl] *= Uav self.dc_energ[icrsh] -= (U_interact + (M-1)*J_hund)/M * 0.5 * Ncr[bl] * Ncr[bl] mpi.report("DC for shell %(icrsh)i and block %(bl)s = %(Uav)f"%locals()) - + # output: mpi.report("DC energy for shell %s = %s"%(icrsh,self.dc_energ[icrsh])) - else: - + else: + a_list = [a for a,al in self.gf_struct_corr[icrsh]] for bl in a_list: self.dc_imp[icrsh][bl] *= use_val - + self.dc_energ[icrsh] = use_val * Ncrtot # output: @@ -743,29 +739,29 @@ class SumkLDA: mpi.report("DC energy = %s"%self.dc_energ[icrsh]) - + def find_dc(self,orb,guess,dens_mat,dens_req=None,precision=0.01): - """searches for DC in order to fulfill charge neutrality. + """Searches for DC in order to fulfill charge neutrality. If dens_req is given, then DC is set such that the LOCAL charge of orbital orb coincides with dens_req.""" - + mu = self.chemical_potential - + def F(dc): self.set_dc(dens_mat=dens_mat,U_interact=0,J_hund=0,orb=orb,use_val=dc) if (dens_req is None): return self.total_density(mu=mu) else: return self.extract_G_loc()[orb].total_density() - + if (dens_req is None): Dens_rel = self.density_required - self.charge_below else: Dens_rel = dens_req - + dcnew = dichotomy.dichotomy(function = F, x_init = guess, y_value = Dens_rel, precision_on_y = precision, delta_x=0.5, @@ -775,16 +771,16 @@ class SumkLDA: return dcnew - - + + def put_Sigma(self, Sigma_imp): """Puts the impurity self energies for inequivalent atoms into the class, respects the multiplicity of the atoms.""" assert isinstance(Sigma_imp,list), "Sigma_imp has to be a list of Sigmas for the correlated shells, even if it is of length 1!" assert len(Sigma_imp)==self.n_inequiv_corr_shells, "give exactly one Sigma for each inequivalent corr. shell!" - + # init self.Sigma_imp: if Sigma_imp[0].note == 'ReFreq': # Real frequency Sigma: @@ -795,7 +791,7 @@ class SumkLDA: # Imaginary frequency Sigma: self.Sigma_imp = [ BlockGf( name_block_generator = [ (a,GfImFreq(indices = al, mesh = Sigma_imp[0].mesh)) for a,al in self.gf_struct_corr[i] ], make_copies = False) for i in xrange(self.n_corr_shells) ] - + # transform the CTQMC blocks to the full matrix: for icrsh in xrange(self.n_corr_shells): s = self.shellmap[icrsh] # s is the index of the inequivalent shell corresponding to icrsh @@ -805,21 +801,19 @@ class SumkLDA: cnt = {} for blname in self.map[s]: cnt[blname] = 0 - - for a,al in self.gf_struct_solver[s]: + + for a,al in self.gf_struct_solver[s].iteritems(): blname = self.map_inv[s][a] map_ind[a] = range(len(al)) for i in al: map_ind[a][i] = cnt[blname] cnt[blname]+=1 - - - for ibl in range(len(self.gf_struct_solver[s])): - for i in range(len(self.gf_struct_solver[s][ibl][1])): - for j in range(len(self.gf_struct_solver[s][ibl][1])): - bl = self.gf_struct_solver[s][ibl][0] - ind1 = self.gf_struct_solver[s][ibl][1][i] - ind2 = self.gf_struct_solver[s][ibl][1][j] + + for bl, orblist in self.gf_struct_solver[s].iteritems(): + for i in range(len(orblist)): + for j in range(len(orblist)): + ind1 = orblist[i] + ind2 = orblist[j] ind1_imp = map_ind[bl][ind1] ind2_imp = map_ind[bl][ind2] self.Sigma_imp[icrsh][self.map_inv[s][bl]][ind1_imp,ind2_imp] <<= Sigma_imp[s][bl][ind1,ind2] @@ -833,17 +827,17 @@ class SumkLDA: def add_dc(self): """Substracts the double counting term from the impurity self energy.""" - + # Be careful: Sigma_imp is already in the global coordinate system!! sres = [s.copy() for s in self.Sigma_imp] for icrsh in xrange(self.n_corr_shells): for bl,gf in sres[icrsh]: dccont = numpy.dot(self.rot_mat[icrsh],numpy.dot(self.dc_imp[icrsh][bl],self.rot_mat[icrsh].conjugate().transpose())) sres[icrsh][bl] -= dccont - - return sres - + return sres + + def set_mu(self,mu): """Sets a new chemical potential""" @@ -893,10 +887,10 @@ class SumkLDA: self.invshellmap = [0] self.n_inequiv_corr_shells = 1 tmp.append( lst[0][1:3] ) - + if (len(lst)>1): for i in range(len(lst)-1): - + fnd = False for j in range(self.n_inequiv_corr_shells): if (tmp[j]==lst[i+1][1:3]): @@ -907,32 +901,32 @@ class SumkLDA: self.n_inequiv_corr_shells += 1 tmp.append( lst[i+1][1:3] ) self.invshellmap.append(i+1) - - + + def total_density(self, mu): """ - Calculates the total charge for the energy window for a given mu. Since in general n_orbitals depends on k, + Calculates the total charge for the energy window for a given mu. Since in general n_orbitals depends on k, the calculation is done in the following order: - G_aa'(k,iw) -> n(k) = Tr G_aa'(k,iw) -> sum_k n_k - + G_aa'(k,iw) -> n(k) = Tr G_aa'(k,iw) -> sum_k n_k + mu: chemical potential - + The calculation is done in the global coordinate system, if distinction is made between local/global! """ - + dens = 0.0 ikarray=numpy.array(range(self.n_k)) - + for ik in mpi.slice_array(ikarray): - - S = self.lattice_gf_matsubara(ik=ik,mu=mu) + + S = self.lattice_gf_matsubara(ik=ik,mu=mu) dens += self.bz_weights[ik] * S.total_density() - + # collect data from mpi: dens = mpi.all_reduce(mpi.world,dens,lambda x,y : x+y) mpi.barrier() - + return dens @@ -944,7 +938,7 @@ class SumkLDA: Dens_rel = self.density_required - self.charge_below - + self.chemical_potential = dichotomy.dichotomy(function = F, x_init = self.chemical_potential, y_value = Dens_rel, precision_on_y = precision, delta_x=0.5, @@ -964,13 +958,13 @@ class SumkLDA: if (orb is None): dens = 0.0 for ish in range(self.n_inequiv_corr_shells): - dens += gnonint[ish].total_density() + dens += gnonint[ish].total_density() else: dens = gnonint[orb].total_density() - + return dens - - + + self.chemical_potential = dichotomy.dichotomy(function = F, x_init = self.chemical_potential, y_value = dens_req, precision_on_y = precision, delta_x=0.5, @@ -982,26 +976,26 @@ class SumkLDA: def extract_G_loc(self, mu=None, with_Sigma = True): - """ + """ extracts the local downfolded Green function at the chemical potential of the class. - At the end, the local G is rotated from the gloabl coordinate system to the local system. + At the end, the local G is rotated from the global coordinate system to the local system. if with_Sigma = False: Sigma is not included => non-interacting local GF """ if (mu is None): mu = self.chemical_potential - - Gloc = [ self.Sigma_imp[icrsh].copy() for icrsh in xrange(self.n_corr_shells) ] # this list will be returned + + Gloc = [ self.Sigma_imp[icrsh].copy() for icrsh in xrange(self.n_corr_shells) ] # this list will be returned for icrsh in xrange(self.n_corr_shells): Gloc[icrsh].zero() # initialize to zero beta = Gloc[0].mesh.beta - + ikarray=numpy.array(range(self.n_k)) - + for ik in mpi.slice_array(ikarray): - - S = self.lattice_gf_matsubara(ik=ik,mu=mu,with_Sigma = with_Sigma, beta = beta) + + S = self.lattice_gf_matsubara(ik=ik,mu=mu,with_Sigma = with_Sigma, beta = beta) S *= self.bz_weights[ik] - + for icrsh in xrange(self.n_corr_shells): tmp = Gloc[icrsh].copy() # init temporary storage for sig,gf in tmp: tmp[sig] <<= self.downfold(ik,icrsh,sig,S[sig],gf) @@ -1012,18 +1006,18 @@ class SumkLDA: Gloc[icrsh] <<= mpi.all_reduce(mpi.world,Gloc[icrsh],lambda x,y : x+y) mpi.barrier() - + # Gloc[:] is now the sum over k projected to the local orbitals. - # here comes the symmetrisation, if needed: + # here comes the symmetrisation, if needed: if (self.symm_op!=0): Gloc = self.Symm_corr.symmetrize(Gloc) - + # Gloc is rotated to the local coordinate system: if (self.use_rotations): for icrsh in xrange(self.n_corr_shells): for sig,gf in Gloc[icrsh]: Gloc[icrsh][sig] <<= self.rotloc(icrsh,gf,direction='toLocal') # transform to CTQMC blocks: - Glocret = [ BlockGf( name_block_generator = [ (a,GfImFreq(indices = al, mesh = Gloc[0].mesh)) for a,al in self.gf_struct_solver[i] ], + Glocret = [ BlockGf( name_block_generator = [ (a,GfImFreq(indices = al, mesh = Gloc[0].mesh)) for a,al in self.gf_struct_solver[i].iteritems() ], make_copies = False) for i in xrange(self.n_inequiv_corr_shells) ] for ish in xrange(self.n_inequiv_corr_shells): @@ -1032,21 +1026,19 @@ class SumkLDA: cnt = {} for blname in self.map[ish]: cnt[blname] = 0 - - for a,al in self.gf_struct_solver[ish]: + + for a,al in self.gf_struct_solver[ish].iteritems(): blname = self.map_inv[ish][a] map_ind[a] = range(len(al)) for i in al: map_ind[a][i] = cnt[blname] cnt[blname]+=1 - - - for ibl in range(len(self.gf_struct_solver[ish])): - for i in range(len(self.gf_struct_solver[ish][ibl][1])): - for j in range(len(self.gf_struct_solver[ish][ibl][1])): - bl = self.gf_struct_solver[ish][ibl][0] - ind1 = self.gf_struct_solver[ish][ibl][1][i] - ind2 = self.gf_struct_solver[ish][ibl][1][j] + + for bl, orblist in self.gf_struct_solver[ish].iteritems(): + for i in range(len(orblist)): + for j in range(len(orblist)): + ind1 = orblist[i] + ind2 = orblist[j] ind1_imp = map_ind[bl][ind1] ind2_imp = map_ind[bl][ind2] Glocret[ish][bl][ind1,ind2] <<= Gloc[self.invshellmap[ish]][self.map_inv[ish][bl]][ind1_imp,ind2_imp] @@ -1054,12 +1046,12 @@ class SumkLDA: # return only the inequivalent shells: return Glocret - + def calc_density_correction(self,filename = 'dens_mat.dat'): """ Calculates the density correction in order to feed it back to the DFT calculations.""" - + assert (type(filename)==StringType), "filename has to be a string!" ntoi = self.names_to_ind[self.SO] @@ -1071,19 +1063,19 @@ class SumkLDA: deltaN[ib] = [ numpy.zeros( [self.n_orbitals[ik,ntoi[ib]],self.n_orbitals[ik,ntoi[ib]]], numpy.complex_) for ik in range(self.n_k)] ikarray=numpy.array(range(self.n_k)) - + dens = {} for ib in bln: dens[ib] = 0.0 - + for ik in mpi.slice_array(ikarray): - + S = self.lattice_gf_matsubara(ik=ik,mu=self.chemical_potential) for sig,g in S: deltaN[sig][ik] = S[sig].density() dens[sig] += self.bz_weights[ik] * S[sig].total_density() - - + + #put mpi Barrier: for sig in deltaN: @@ -1092,7 +1084,7 @@ class SumkLDA: dens[sig] = mpi.all_reduce(mpi.world,dens[sig],lambda x,y : x+y) mpi.barrier() - + # now save to file: if (mpi.is_master_node()): if (self.SP==0): @@ -1151,6 +1143,6 @@ class SumkLDA: f1.write("\n") f1.write("\n") f1.close() - + return deltaN, dens diff --git a/python/sumk_lda_tools.py b/python/sumk_lda_tools.py index 1f3faa53..e7ee4947 100644 --- a/python/sumk_lda_tools.py +++ b/python/sumk_lda_tools.py @@ -24,23 +24,17 @@ from types import * import numpy import pytriqs.utility.dichotomy as dichotomy from pytriqs.gf.local import * -#from pytriqs.applications.impurity_solvers.operators import * from pytriqs.operators import * import pytriqs.utility.mpi as mpi from datetime import datetime - -#from pytriqs.applications.dft.symmetry import * -#from pytriqs.applications.dft.sumk_lda import SumkLDA from symmetry import * from sumk_lda import SumkLDA - import string - def read_fortran_file (filename): """ Returns a generator that yields all numbers in the Fortran file as float, one by one""" import os.path - if not(os.path.exists(filename)) : raise IOError, "File %s does not exists"%filename + if not(os.path.exists(filename)) : raise IOError, "File %s does not exist."%filename for line in open(filename,'r') : for x in line.replace('D','E').split() : yield string.atof(x) @@ -56,12 +50,12 @@ class SumkLDATools(SumkLDA): self.Gupf_refreq = None SumkLDA.__init__(self,hdf_file=hdf_file,mu=mu,h_field=h_field,use_lda_blocks=use_lda_blocks,lda_data=lda_data, symm_corr_data=symm_corr_data,par_proj_data=par_proj_data,symm_par_data=symm_par_data, - bands_data=bands_data) - + bands_data=bands_data) + def downfold_pc(self,ik,ir,ish,sig,gf_to_downfold,gf_inp): """Downfolding a block of the Greens function""" - + gf_downfolded = gf_inp.copy() isp = self.names_to_ind[self.SO][sig] # get spin index for proj. matrices dim = self.shells[ish][3] @@ -87,15 +81,15 @@ class SumkLDATools(SumkLDA): gf_rotated.from_L_G_R(self.rot_mat_all[ish].conjugate(),gf_rotated,self.rot_mat_all[ish].transpose()) else: gf_rotated.from_L_G_R(self.rot_mat_all[ish],gf_rotated,self.rot_mat_all[ish].conjugate().transpose()) - + elif (direction=='toLocal'): if ((self.rot_mat_all_time_inv[ish]==1)and(self.SO)): gf_rotated <<= gf_rotated.transpose() gf_rotated.from_L_G_R(self.rot_mat_all[ish].transpose(),gf_rotated,self.rot_mat_all[ish].conjugate()) else: gf_rotated.from_L_G_R(self.rot_mat_all[ish].conjugate().transpose(),gf_rotated,self.rot_mat_all[ish]) - - + + return gf_rotated @@ -107,7 +101,7 @@ class SumkLDATools(SumkLDA): bln = self.block_names[self.SO] if (not hasattr(self,"Sigma_imp")): with_Sigma=False - if (with_Sigma): + if (with_Sigma): assert self.Sigma_imp[0].note == 'ReFreq', "Real frequency Sigma needed for lattice_gf_realfreq!" stmp = self.add_dc() else: @@ -139,7 +133,7 @@ class SumkLDATools(SumkLDA): glist = lambda : [ GfReFreq(indices = al, window=(mesh[0],mesh[1]),n_points=mesh[2]) for a,al in gf_struct] self.Gupf_refreq = BlockGf(name_list = a_list, block_list = glist(),make_copies=False) self.Gupf_refreq.zero() - + idmat = [numpy.identity(self.n_orbitals[ik,ntoi[bl]],numpy.complex_) for bl in bln] self.Gupf_refreq <<= Omega + 1j*broadening @@ -163,7 +157,7 @@ class SumkLDATools(SumkLDA): def check_input_dos(self, om_min, om_max, n_om, beta=10, broadening=0.01): - + delta_om = (om_max-om_min)/(n_om-1) om_mesh = numpy.zeros([n_om],numpy.float_) @@ -189,59 +183,59 @@ class SumkLDATools(SumkLDA): glist = lambda : [ GfReFreq(indices = al, window = (om_min,om_max), n_points = n_om) for a,al in self.gf_struct_corr[icrsh]] Gloc.append(BlockGf(name_list = b_list, block_list = glist(),make_copies=False)) for icrsh in xrange(self.n_corr_shells): Gloc[icrsh].zero() # initialize to zero - + for ik in xrange(self.n_k): Gupf=self.lattice_gf_realfreq(ik=ik,mu=self.chemical_potential,broadening=broadening,mesh=(om_min,om_max,n_om),with_Sigma=False) Gupf *= self.bz_weights[ik] # non-projected DOS - for iom in range(n_om): - for sig,gf in Gupf: + for iom in range(n_om): + for sig,gf in Gupf: asd = gf.data[iom,:,:].imag.trace()/(-3.1415926535) DOS[sig][iom] += asd - + for icrsh in xrange(self.n_corr_shells): tmp = Gloc[icrsh].copy() for sig,gf in tmp: tmp[sig] <<= self.downfold(ik,icrsh,sig,Gupf[sig],gf) # downfolding G Gloc[icrsh] += tmp - - + + if (self.symm_op!=0): Gloc = self.Symm_corr.symmetrize(Gloc) if (self.use_rotations): for icrsh in xrange(self.n_corr_shells): for sig,gf in Gloc[icrsh]: Gloc[icrsh][sig] <<= self.rotloc(icrsh,gf,direction='toLocal') - + # Gloc can now also be used to look at orbitally resolved quantities for ish in range(self.n_inequiv_corr_shells): for sig,gf in Gloc[self.invshellmap[ish]]: # loop over spins - for iom in range(n_om): DOSproj[ish][sig][iom] += gf.data[iom,:,:].imag.trace()/(-3.1415926535) + for iom in range(n_om): DOSproj[ish][sig][iom] += gf.data[iom,:,:].imag.trace()/(-3.1415926535) DOSproj_orb[ish][sig][:,:,:] += gf.data[:,:,:].imag/(-3.1415926535) - + # output: if (mpi.is_master_node()): for bn in self.block_names[self.SO]: f=open('DOS%s.dat'%bn, 'w') for i in range(n_om): f.write("%s %s\n"%(om_mesh[i],DOS[bn][i])) - f.close() + f.close() for ish in range(self.n_inequiv_corr_shells): f=open('DOS%s_proj%s.dat'%(bn,ish),'w') for i in range(n_om): f.write("%s %s\n"%(om_mesh[i],DOSproj[ish][bn][i])) - f.close() - + f.close() + for i in range(self.corr_shells[self.invshellmap[ish]][3]): for j in range(i,self.corr_shells[self.invshellmap[ish]][3]): Fname = 'DOS'+bn+'_proj'+str(ish)+'_'+str(i)+'_'+str(j)+'.dat' f=open(Fname,'w') for iom in range(n_om): f.write("%s %s\n"%(om_mesh[iom],DOSproj_orb[ish][bn][iom,i,j])) f.close() - - + + def read_par_proj_input_from_hdf(self): """ @@ -268,7 +262,7 @@ class SumkLDATools(SumkLDA): mu = self.chemical_potential gf_struct_proj = [ [ (al, range(self.shells[i][3])) for al in self.block_names[self.SO] ] for i in xrange(self.n_shells) ] - Gproj = [BlockGf(name_block_generator = [ (a,GfReFreq(indices = al, mesh = self.Sigma_imp[0].mesh)) for a,al in gf_struct_proj[ish] ], make_copies = False ) + Gproj = [BlockGf(name_block_generator = [ (a,GfReFreq(indices = al, mesh = self.Sigma_imp[0].mesh)) for a,al in gf_struct_proj[ish] ], make_copies = False ) for ish in xrange(self.n_shells)] for ish in range(self.n_shells): Gproj[ish].zero() @@ -295,32 +289,32 @@ class SumkLDATools(SumkLDA): S *= self.bz_weights[ik] # non-projected DOS - for iom in range(n_om): + for iom in range(n_om): for sig,gf in S: DOS[sig][iom] += gf.data[iom,:,:].imag.trace()/(-3.1415926535) - + #projected DOS: for ish in xrange(self.n_shells): tmp = Gproj[ish].copy() for ir in xrange(self.n_parproj[ish]): for sig,gf in tmp: tmp[sig] <<= self.downfold_pc(ik,ir,ish,sig,S[sig],gf) Gproj[ish] += tmp - + # collect data from mpi: for sig in DOS: DOS[sig] = mpi.all_reduce(mpi.world,DOS[sig],lambda x,y : x+y) for ish in xrange(self.n_shells): Gproj[ish] <<= mpi.all_reduce(mpi.world,Gproj[ish],lambda x,y : x+y) - mpi.barrier() - + mpi.barrier() + if (self.symm_op!=0): Gproj = self.Symm_par.symmetrize(Gproj) # rotation to local coord. system: if (self.use_rotations): for ish in xrange(self.n_shells): for sig,gf in Gproj[ish]: Gproj[ish][sig] <<= self.rotloc_all(ish,gf,direction='toLocal') - + for ish in range(self.n_shells): - for sig,gf in Gproj[ish]: + for sig,gf in Gproj[ish]: for iom in range(n_om): DOSproj[ish][sig][iom] += gf.data[iom,:,:].imag.trace()/(-3.1415926535) DOSproj_orb[ish][sig][:,:,:] += gf.data[:,:,:].imag / (-3.1415926535) @@ -330,14 +324,14 @@ class SumkLDATools(SumkLDA): for bn in self.block_names[self.SO]: f=open('./DOScorr%s.dat'%bn, 'w') for i in range(n_om): f.write("%s %s\n"%(Msh[i],DOS[bn][i])) - f.close() + f.close() # partial for ish in range(self.n_shells): f=open('DOScorr%s_proj%s.dat'%(bn,ish),'w') for i in range(n_om): f.write("%s %s\n"%(Msh[i],DOSproj[ish][bn][i])) f.close() - + for i in range(self.shells[ish][3]): for j in range(i,self.shells[ish][3]): Fname = './DOScorr'+bn+'_proj'+str(ish)+'_'+str(i)+'_'+str(j)+'.dat' @@ -349,7 +343,7 @@ class SumkLDATools(SumkLDA): def spaghettis(self,broadening,shift=0.0,plot_range=None, ishell=None, invert_Akw=False, fermi_surface=False): - """ Calculates the correlated band structure with a real-frequency self energy. + """ Calculates the correlated band structure with a real-frequency self energy. ATTENTION: Many things from the original input file are are overwritten!!!""" assert hasattr(self,"Sigma_imp"), "Set Sigma First!!" @@ -358,13 +352,13 @@ class SumkLDATools(SumkLDA): if not retval: return retval if fermi_surface: ishell=None - + # print hamiltonian for checks: if ((self.SP==1)and(self.SO==0)): f1=open('hamup.dat','w') f2=open('hamdn.dat','w') - - for ik in xrange(self.n_k): + + for ik in xrange(self.n_k): for i in xrange(self.n_orbitals[ik,0]): f1.write('%s %s\n'%(ik,self.hopping[ik,0,i,i].real)) for i in xrange(self.n_orbitals[ik,1]): @@ -381,7 +375,7 @@ class SumkLDATools(SumkLDA): f.write('\n') f.close() - + #========================================= # calculate A(k,w): @@ -419,17 +413,17 @@ class SumkLDATools(SumkLDA): for ik in xrange(self.n_k): - S = self.lattice_gf_realfreq(ik=ik,mu=mu,broadening=broadening) + S = self.lattice_gf_realfreq(ik=ik,mu=mu,broadening=broadening) if (ishell is None): # non-projected A(k,w) - for iom in range(n_om): + for iom in range(n_om): if (M[iom]>om_minplot) and (M[iom]om_minplot) and (M[iom]om_minplot) and (M[iom]om_minplot) and (M[iom]