diff --git a/python/block_structure.py b/python/block_structure.py index b0903161..f1501409 100644 --- a/python/block_structure.py +++ b/python/block_structure.py @@ -2,12 +2,13 @@ import copy import numpy as np from pytriqs.gf import GfImFreq, BlockGf from ast import literal_eval +import pytriqs.utility.mpi as mpi from warnings import warn class BlockStructure(object): """ Contains information about the Green function structure. - This class contains information about the structure of the solver + This class contains information about the structure of the solver and sumk Green functions and the mapping between them. Parameters @@ -33,19 +34,21 @@ class BlockStructure(object): solver_to_sumk_block : list of dict solver_to_sumk_block[ish][from_block] = to_block - maps from the solver block to the sumk block + maps from the solver block to the sumk block for *inequivalent* correlated shell ish """ def __init__(self,gf_struct_sumk=None, gf_struct_solver=None, solver_to_sumk=None, sumk_to_solver=None, - solver_to_sumk_block=None): + solver_to_sumk_block=None, + deg_shells=None): self.gf_struct_sumk = gf_struct_sumk self.gf_struct_solver = gf_struct_solver self.solver_to_sumk = solver_to_sumk self.sumk_to_solver = sumk_to_solver self.solver_to_sumk_block = solver_to_sumk_block + self.deg_shells = deg_shells @classmethod def full_structure(cls,gf_struct,corr_to_inequiv): @@ -99,20 +102,21 @@ class BlockStructure(object): gf_struct_sumk = gs_sumk_all, solver_to_sumk = copy.deepcopy(solver_to_sumk), sumk_to_solver = solver_to_sumk, - solver_to_sumk_block = s2sblock) + solver_to_sumk_block = s2sblock, + deg_shells = [[] for ish in range(len(gf_struct))]) def pick_gf_struct_solver(self,new_gf_struct): - """ Pick selected orbitals within blocks. + """ Pick selected orbitals within blocks. This throws away parts of the Green's function that (for some - reason - be sure that you know what you're doing) shouldn't be + reason - be sure that you know what you're doing) shouldn't be included in the calculation. To drop an entire block, just don't include it. To drop a certain index within a block, just don't include it. - If it was before: - + If it was before: + [{'up':[0,1],'down':[0,1],'left':[0,1]}] to choose the 0th index of the up block and the 1st index of @@ -130,11 +134,11 @@ class BlockStructure(object): Parameters ---------- new_gf_struct : list of dict - formatted the same as gf_struct_solver: + formatted the same as gf_struct_solver: new_gf_struct[ish][block]=list of indices in that block. """ - + for ish in range(len(self.gf_struct_solver)): gf_struct = new_gf_struct[ish] @@ -154,24 +158,24 @@ class BlockStructure(object): new_ind = gf_struct[blk].index(ind) self.sumk_to_solver[ish][k]=(blk,new_ind) else: - self.sumk_to_solver[ish][k]=(None,None) + self.sumk_to_solver[ish][k]=(None,None) # reindexing gf_struct so that it starts with 0 for k in gf_struct: gf_struct[k]=range(len(gf_struct[k])) self.gf_struct_solver[ish]=gf_struct def pick_gf_struct_sumk(self,new_gf_struct): - """ Pick selected orbitals within blocks. + """ Pick selected orbitals within blocks. This throws away parts of the Green's function that (for some - reason - be sure that you know what you're doing) shouldn't be + reason - be sure that you know what you're doing) shouldn't be included in the calculation. To drop an entire block, just don't include it. To drop a certain index within a block, just don't include it. - If it was before: - + If it was before: + [{'up':[0,1],'down':[0,1],'left':[0,1]}] to choose the 0th index of the up block and the 1st index of @@ -188,11 +192,11 @@ class BlockStructure(object): Parameters ---------- new_gf_struct : list of dict - formatted the same as gf_struct_solver: + formatted the same as gf_struct_solver: new_gf_struct[ish][block]=list of indices in that block. - However, the indices are not according to the solver Gf + However, the indices are not according to the solver Gf but the sumk Gf. """ @@ -218,7 +222,7 @@ class BlockStructure(object): Parameters ---------- mapping : list of dict - the dict consists of elements + the dict consists of elements (from_block,from_index) : (to_block,to_index) that maps from one structure to the other """ @@ -254,7 +258,7 @@ class BlockStructure(object): def create_gf(self,ish=0,gf_function=GfImFreq,**kwargs): """ Create a zero BlockGf having the gf_struct_solver structure. - When using GfImFreq as gf_function, typically you have to + When using GfImFreq as gf_function, typically you have to supply beta as keyword argument. Parameters @@ -284,7 +288,7 @@ class BlockStructure(object): .. warning:: Elements that are zero in the new structure due to - the new block structure will be just ignored, thus + the new block structure will be just ignored, thus approximated to zero. Parameters @@ -292,15 +296,24 @@ class BlockStructure(object): G : BlockGf the Gf that should be converted G_struct : GfStructure - the structure ofthat G + the structure of that G ish : int shell index - show_warnings : bool - whether to show warnings when elements of the Green's + show_warnings : bool or float + whether to show warnings when elements of the Green's function get thrown away + if float, set the threshold for the magnitude of an element + about to be thrown away to trigger a warning + (default: 1.e-10) **kwargs : options passed to the constructor for the new Gf """ + + warning_threshold = 1.e-10 + if isinstance(show_warnings, float): + warning_threshold = show_warnings + show_warnings = True + G_new = self.create_gf(ish=ish,**kwargs) for block in G_struct.gf_struct_solver[ish].keys(): for i1 in G_struct.gf_struct_solver[ish][block]: @@ -311,22 +324,24 @@ class BlockStructure(object): i2_sol = self.sumk_to_solver[ish][i2_sumk] if i1_sol[0] is None or i2_sol[0] is None: if show_warnings: - warn(('Element {},{} of block {} of G is not present '+ - 'in the new structure').format(i1,i2,block)) + if mpi.is_master_node(): + warn(('Element {},{} of block {} of G is not present '+ + 'in the new structure').format(i1,i2,block)) continue if i1_sol[0]!=i2_sol[0]: - if show_warnings: - warn(('Element {},{} of block {} of G is approximated '+ - 'to zero to match the new structure.').format( - i1,i2,block)) + if show_warnings and np.max(np.abs(G[block][i1,i2].data)) > warning_threshold: + if mpi.is_master_node(): + warn(('Element {},{} of block {} of G is approximated '+ + 'to zero to match the new structure. Max abs value: {}').format( + i1,i2,block,np.max(np.abs(G[block][i1,i2].data)))) continue G_new[i1_sol[0]][i1_sol[1],i2_sol[1]] = \ G[block][i1,i2] return G_new def approximate_as_diagonal(self): - """ Create a structure for a GF with zero off-diagonal elements. - + """ Create a structure for a GF with zero off-diagonal elements. + .. warning:: In general, this will throw away non-zero elements of the @@ -351,7 +366,8 @@ class BlockStructure(object): def __eq__(self,other): def compare(one,two): if type(one)!=type(two): - return False + if not (isinstance(one, (bool, np.bool_)) and isinstance(two, (bool, np.bool_))): + return False if one is None and two is None: return True if isinstance(one,list) or isinstance(one,tuple): @@ -361,10 +377,10 @@ class BlockStructure(object): if not compare(x,y): return False return True - elif isinstance(one,int): - return one==two - elif isinstance(one,str): + elif isinstance(one,(int,bool, str, np.bool_)): return one==two + elif isinstance(one,np.ndarray): + return np.all(one==two) elif isinstance(one,dict): if set(one.keys()) != set(two.keys()): return False @@ -375,8 +391,9 @@ class BlockStructure(object): warn('Cannot compare {}'.format(type(one))) return False - for prop in [ "gf_struct_sumk", "gf_struct_solver", - "solver_to_sumk", "sumk_to_solver", "solver_to_sumk_block"]: + for prop in [ "gf_struct_sumk", "gf_struct_solver", + "solver_to_sumk", "sumk_to_solver", "solver_to_sumk_block", + "deg_shells"]: if not compare(getattr(self,prop),getattr(other,prop)): return False return True @@ -388,8 +405,8 @@ class BlockStructure(object): """ Reduce to dict for HDF5 export.""" ret = {} - for element in [ "gf_struct_sumk", "gf_struct_solver", - "solver_to_sumk_block"]: + for element in [ "gf_struct_sumk", "gf_struct_solver", + "solver_to_sumk_block","deg_shells"]: ret[element] = getattr(self,element) def construct_mapping(mapping): @@ -436,6 +453,18 @@ class BlockStructure(object): keys = sorted(element[ish].keys(),key=keyfun) for k in keys: s+=' '+str(k)+str(element[ish][k])+'\n' + s += "deg_shells\n" + for ish in range(len(self.deg_shells)): + s+=' shell '+str(ish)+'\n' + for l in range(len(self.deg_shells[ish])): + s+=' equivalent group '+str(l)+'\n' + if isinstance(self.deg_shells[ish][l],dict): + for key, val in self.deg_shells[ish][l].iteritems(): + s+=' '+key+('*' if val[1] else '')+':\n' + s+=' '+str(val[0]).replace('\n','\n ')+'\n' + else: + for key in self.deg_shells[ish][l]: + s+=' '+key+'\n' return s from pytriqs.archive.hdf_archive_schemes import register_class diff --git a/python/sumk_dft.py b/python/sumk_dft.py index 900faff5..8f56f544 100644 --- a/python/sumk_dft.py +++ b/python/sumk_dft.py @@ -25,12 +25,15 @@ import numpy import pytriqs.utility.dichotomy as dichotomy from pytriqs.gf import * import pytriqs.utility.mpi as mpi +from pytriqs.utility.comparison_tests import assert_arrays_are_close from pytriqs.archive import * from symmetry import * from block_structure import BlockStructure from sets import Set from itertools import product from warnings import warn +from scipy import compress +from scipy.optimize import minimize class SumkDFT(object): @@ -848,6 +851,414 @@ class SumkDFT(object): elif (ind1 < 0) and (ind2 < 0): self.deg_shells[ish].append([block1, block2]) + def _get_hermitian_quantity_from_gf(self, G): + """ Convert G to a Hermitian quantity + + For G(tau) and G(iw), G(tau) is returned. + For G(t) and G(w), the spectral function is returned. + + Parameters + ---------- + G : list of BlockGf of GfImFreq, GfImTime, GfReFreq or GfReTime + the input Green's function + + Returns + ------- + gf : list of BlockGf of GfImTime or GfReFreq + the output G(tau) or A(w) + """ + # make a GfImTime from the supplied GfImFreq + if all(isinstance(g_sh._first(), GfImFreq) for g_sh in G): + gf = [BlockGf(name_block_generator = [(name, GfImTime(beta=block.mesh.beta, + indices=block.indices,n_points=len(block.mesh)+1)) for name, block in g_sh]) + for g_sh in G] + for ish in range(len(gf)): + for name, g in gf[ish]: + g.set_from_inverse_fourier(G[ish][name]) + # keep a GfImTime from the supplied GfImTime + elif all(isinstance(g_sh._first(), GfImTime) for g_sh in G): + gf = G + # make a spectral function from the supplied GfReFreq + elif all(isinstance(g_sh._first(), GfReFreq) for g_sh in G): + gf = [g_sh.copy() for g_sh in G] + for ish in range(len(gf)): + for name, g in gf[ish]: + g << 1.0j*(g-g.conjugate().transpose())/2.0/numpy.pi + elif all(isinstance(g_sh._first(), GfReTime) for g_sh in G): + def get_delta_from_mesh(mesh): + w0 = None + for w in mesh: + if w0 is None: + w0 = w + else: + return w-w0 + gf = [BlockGf( + name_block_generator = [(name, + GfReFreq( + window=(-numpy.pi*(len(block.mesh)-1) / (len(block.mesh)*get_delta_from_mesh(block.mesh)), numpy.pi*(len(block.mesh)-1) / (len(block.mesh)*get_delta_from_mesh(block.mesh))), + n_points=len(block.mesh), + indices=block.indices)) for name, block in g_sh]) + for g_sh in G] + + for ish in range(len(gf)): + for name, g in gf[ish]: + g.set_from_fourier(G[ish][name]) + g << 1.0j*(g-g.conjugate().transpose())/2.0/numpy.pi + else: + raise Exception("G must be a list of BlockGf of either GfImFreq, GfImTime, GfReFreq or GfReTime") + return gf + + + + def analyse_block_structure_from_gf(self, G, threshold=1.e-5, include_shells=None, analyse_deg_shells = True): + r""" + Determines the block structure of local Green's functions by analysing + the structure of the corresponding non-interacting Green's function. + The resulting block structures for correlated shells are + stored in the :class:`SumkDFT.block_structure ` + attribute. + + This is a safer alternative to analyse_block_structure, because + the full non-interacting Green's function is taken into account + and not just the density matrix and Hloc. + + Parameters + ---------- + G : list of BlockGf of GfImFreq, GfImTime, GfReFreq or GfReTime + the non-interacting Green's function for each inequivalent correlated shell + threshold : real, optional + If the difference between matrix elements is below threshold, + they are considered to be equal. + include_shells : list of integers, optional + List of correlated shells to be analysed. + If include_shells is not provided all correlated shells will be analysed. + analyse_deg_shells : bool + Whether to call the analyse_deg_shells function + after having finished the block structure analysis + + Returns + ------- + G : list of BlockGf of GfImFreq or GfImTime + the Green's function transformed into the new block structure + """ + + gf = self._get_hermitian_quantity_from_gf(G) + + # initialize the variables + self.gf_struct_solver = [{} for ish in range(self.n_inequiv_shells)] + self.sumk_to_solver = [{} for ish in range(self.n_inequiv_shells)] + self.solver_to_sumk = [{} for ish in range(self.n_inequiv_shells)] + self.solver_to_sumk_block = [{} + for ish in range(self.n_inequiv_shells)] + + # the maximum value of each matrix element of each block and shell + max_gf = [{name:numpy.max(numpy.abs(g.data),0) for name, g in gf[ish]} for ish in range(self.n_inequiv_shells)] + + if include_shells is None: + # include all shells + include_shells = range(self.n_inequiv_shells) + + for ish in include_shells: + for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]]['SO']]: + n_orb = self.corr_shells[self.inequiv_to_corr[ish]]['dim'] + + # gives an index list of entries larger that threshold + maxgf_bool = (abs(max_gf[ish][sp]) > threshold) + + # Determine off-diagonal entries in upper triangular part of the + # Green's function + offdiag = Set([]) + for i in range(n_orb): + for j in range(i + 1, n_orb): + if maxgf_bool[i, j]: + offdiag.add((i, j)) + + # Determine the number of non-hybridising blocks in the gf + blocs = [[i] for i in range(n_orb)] + while len(offdiag) != 0: + pair = offdiag.pop() + for b1, b2 in product(blocs, blocs): + if (pair[0] in b1) and (pair[1] in b2): + if blocs.index(b1) != blocs.index(b2): # In separate blocks? + # Merge two blocks + b1.extend(blocs.pop(blocs.index(b2))) + break # Move on to next pair in offdiag + + # Set the gf_struct for the solver accordingly + num_blocs = len(blocs) + for i in range(num_blocs): + blocs[i].sort() + self.gf_struct_solver[ish].update( + [('%s_%s' % (sp, i), range(len(blocs[i])))]) + + # Construct sumk_to_solver taking (sumk_block, sumk_index) --> (solver_block, solver_inner) + # and solver_to_sumk taking (solver_block, solver_inner) --> + # (sumk_block, sumk_index) + for i in range(num_blocs): + for j in range(len(blocs[i])): + block_sumk = sp + inner_sumk = blocs[i][j] + block_solv = '%s_%s' % (sp, i) + inner_solv = j + self.sumk_to_solver[ish][(block_sumk, inner_sumk)] = ( + block_solv, inner_solv) + self.solver_to_sumk[ish][(block_solv, inner_solv)] = ( + block_sumk, inner_sumk) + self.solver_to_sumk_block[ish][block_solv] = block_sumk + + # transform G to the new structure + full_structure = BlockStructure.full_structure( + [{sp:range(self.corr_shells[self.inequiv_to_corr[ish]]['dim']) + for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]]['SO']]} + for ish in range(self.n_inequiv_shells)],None) + G_transformed = [ + self.block_structure.convert_gf(G[ish], + full_structure, ish, mesh=G[ish].mesh.copy(), show_warnings=threshold, + gf_function=type(G[ish]._first())) + for ish in range(self.n_inequiv_shells)] + + if analyse_deg_shells: + self.analyse_deg_shells(G_transformed, threshold, include_shells) + return G_transformed + + def analyse_deg_shells(self, G, threshold=1.e-5, include_shells=None): + r""" + Determines the degenerate shells of local Green's functions by analysing + the structure of the corresponding non-interacting Green's function. + The results are stored in the + :class:`SumkDFT.block_structure ` + attribute. + + Due to the implementation and numerics, the maximum difference between + two matrix elements that are detected as equal can be a bit higher + (e.g. a factor of two) than the actual threshold. + + Parameters + ---------- + G : list of BlockGf of GfImFreq or GfImTime + the non-interacting Green's function for each inequivalent correlated shell + threshold : real, optional + If the difference between matrix elements is below threshold, + they are considered to be equal. + include_shells : list of integers, optional + List of correlated shells to be analysed. + If include_shells is not provided all correlated shells will be analysed. + """ + + # initialize + self.deg_shells = [[] for ish in range(self.n_inequiv_shells)] + + # helper function + def null(A, eps=1e-15): + """ Calculate the null-space of matrix A """ + u, s, vh = numpy.linalg.svd(A) + null_mask = (s <= eps) + null_space = compress(null_mask, vh, axis=0) + return null_space.conjugate().transpose() + + gf = self._get_hermitian_quantity_from_gf(G) + + if include_shells is None: + # include all shells + include_shells = range(self.n_inequiv_shells) + + # We consider two blocks equal, if their Green's functions obey + # maybe_conjugate1( v1^dagger G1 v1 ) = maybe_conjugate2( v2^dagger G2 v2 ) + # where maybe_conjugate is a function that conjugates the Green's + # function if the flag 'conjugate' is set and the v are unitary + # matrices + # + # for each pair of blocks, we check whether there is a transformation + # maybe_conjugate( T G1 T^dagger ) = G2 + # where our goal is to find T + # we just try whether there is such a T with and without conjugation + for ish in include_shells: + for block1 in self.gf_struct_solver[ish].iterkeys(): + for block2 in self.gf_struct_solver[ish].iterkeys(): + if block1==block2: continue + + # check if the blocks are already present in the deg_shells + ind1 = -1 + ind2 = -2 + for n, ind in enumerate(self.deg_shells[ish]): + if block1 in ind: + ind1 = n + v1 = ind[block1] + if block2 in ind: + ind2 = n + v2 = ind[block2] + + # if both are already present, go to the next pair of blocks + if ind1 >= 0 and ind2 >= 0: + continue + + gf1 = gf[ish][block1] + gf2 = gf[ish][block2] + + # the two blocks have to have the same shape + if gf1.target_shape != gf2.target_shape: + continue + + # Instead of directly comparing the two blocks, we + # compare its eigenvalues. As G(tau) is Hermitian, + # they are real and the eigenvector matrix is unitary. + # Thus, if the eigenvalues are equal we can transform + # one block to make it equal to the other (at least + # for tau=0). + + e1 = numpy.linalg.eigvalsh(gf1.data[0]) + e2 = numpy.linalg.eigvalsh(gf2.data[0]) + if numpy.any(abs(e1-e2) > threshold): continue + + for conjugate in [False,True]: + if conjugate: + gf2 = gf2.conjugate() + + # we want T gf1 T^dagger = gf2 + # while for a given tau, T could be calculated + # by diagonalizing gf1 and gf2, this does not + # work for all taus simultaneously because of + # numerical imprecisions + + # rather, we rewrite the equation to + # T gf1 = gf2 T + # which is the Sylvester equation. + # For that equation, one can use the Kronecker + # product to get a linear problem, which consists + # of finding the null space of M vec T = 0. + + M = numpy.kron(numpy.eye(*gf1.target_shape),gf2.data[0])-numpy.kron(gf1.data[0].transpose(),numpy.eye(*gf1.target_shape)) + N = null(M, threshold) + + # now we get the intersection of the null spaces + # of all values of tau + for i in range(1,len(gf1.data)): + M = numpy.kron(numpy.eye(*gf1.target_shape),gf2.data[i])-numpy.kron(gf1.data[i].transpose(),numpy.eye(*gf1.target_shape)) + # transform M into current null space + M = numpy.dot(M, N) + N = numpy.dot(N, null(M, threshold)) + if numpy.size(N) == 0: + break + + # no intersection of the null spaces -> no symmetry + if numpy.size(N) == 0: continue + + # reshape N: it then has the indices matrix, matrix, number of basis vectors of the null space + N = N.reshape(gf1.target_shape[0], gf1.target_shape[1], -1).transpose([1, 0, 2]) + + """ + any matrix in the null space can now be constructed as + M = 0 + for i in range(N.shape[-1]): + M += y[i]*N[:,:,i] + with coefficients (complex numbers) y[i]. + + We want to get a set of coefficients y so that M is unitary. + Unitary means M M^dagger = 1. + Thus, + sum y[i] N[:,:,i] y[j].conjugate() N[:,:,j].conjugate().transpose() = eye. + The object N[:,:,i] N[:,:,j] is a four-index object which we call Z. + """ + Z = numpy.einsum('aci,bcj->abij', N, N.conjugate()) + + """ + function chi2 + This function takes a real parameter vector y and reinterprets it as complex. + Then, it calculates the chi2 of + sum y[i] N[:,:,i] y[j].conjugate() N[:,:,j].conjugate().transpose() - eye. + """ + def chi2(y): + # reinterpret y as complex number + y = y.view(numpy.complex_) + ret = 0.0 + for a in range(Z.shape[0]): + for b in range(Z.shape[1]): + ret += numpy.abs(numpy.dot(y, numpy.dot(Z[a, b], y.conjugate())) + - (1.0 if a == b else 0.0))**2 + return ret + + # use the minimization routine from scipy + res = minimize(chi2, numpy.ones(2 * N.shape[-1])) + + # if the minimization fails, there is probably no symmetry + if not res.success: continue + # check if the minimization returned zero within the tolerance + if res.fun > threshold: continue + + # reinterpret the solution as a complex number + y = res.x.view(numpy.complex_) + + # reconstruct the T matrix + T = numpy.zeros(N.shape[:-1], dtype=numpy.complex_) + for i in range(len(y)): + T += N[:, :, i] * y[i] + + # transform gf1 using T + G_transformed = gf1.copy() + G_transformed.from_L_G_R(T, gf1, T.conjugate().transpose()) + + # it does not make sense to check the tails for an + # absolute error because it will usually not hold; + # we could just check the relative error + # (here, we ignore it, reasoning that if the data + # is the same, the tails have to coincide as well) + try: + assert_arrays_are_close(G_transformed.data, gf2.data, threshold) + except (RuntimeError, AssertionError): + # the symmetry does not hold + continue + + # Now that we have found a valid T, we have to + # rewrite it to match the convention that + # C1(v1^dagger G1 v1) = C2(v2^dagger G2 v2), + # where C conjugates if the flag is True + + # For each group of degenerate shells, the list + # SK.deg_shells[ish] contains a dict. The keys + # of the dict are the block names, the values + # are tuples. The first entry of the tuple is + # the transformation matrix v, the second entry + # is the conjugation flag + + # the second block is already present + # set v1 and C1 so that they are compatible with + # C(T gf1 T^dagger) = gf2 + # and with + # C1(v1^dagger G1 v1) = C2(v2^dagger G2 v2) + if (ind1 < 0) and (ind2 >= 0): + if conjugate: + self.deg_shells[ish][ind2][block1] = numpy.dot(T.conjugate().transpose(), v2[0].conjugate()), not v2[1] + else: + self.deg_shells[ish][ind2][block1] = numpy.dot(T.conjugate().transpose(), v2[0]), v2[1] + # the first block is already present + # set v2 and C2 so that they are compatible with + # C(T gf1 T^dagger) = gf2 + # and with + # C1(v1^dagger G1 v1) = C2(v2^dagger G2 v2) + elif (ind1 >= 0) and (ind2 < 0): + if conjugate: + self.deg_shells[ish][ind1][block2] = numpy.dot(T.conjugate(), v1[0].conjugate()), not v1[1] + else: + self.deg_shells[ish][ind1][block2] = numpy.dot(T, v1[0]), v1[1] + # the blocks are not already present + # we arbitrarily choose v1=eye and C1=False and + # set v2 and C2 so that they are compatible with + # C(T gf1 T^dagger) = gf2 + # and with + # C1(v1^dagger G1 v1) = C2(v2^dagger G2 v2) + elif (ind1 < 0) and (ind2 < 0): + d = dict() + d[block1] = numpy.eye(*gf1.target_shape), False + if conjugate: + d[block2] = T.conjugate(), True + else: + d[block2] = T, False + self.deg_shells[ish].append(d) + + # a block was found, break out of the loop + break + + def density_matrix(self, method='using_gf', beta=40.0): """Calculate density matrices in one of two ways. @@ -1212,20 +1623,52 @@ class SumkDFT(object): Parameters ---------- gf_to_symm : gf_struct_solver like - Input GF. + Input and output GF (i.e., it gets overwritten) orb : int Index of an inequivalent shell. """ + # when reading block_structures written with older versions from + # an h5 file, self.deg_shells might be None + if self.deg_shells is None: return + for degsh in self.deg_shells[orb]: - ss = gf_to_symm[degsh[0]].copy() - ss.zero() + # ss will hold the averaged orbitals in the basis where the + # blocks are all equal + # i.e. maybe_conjugate(v^dagger gf v) + ss = None n_deg = len(degsh) - for bl in degsh: - ss += gf_to_symm[bl] / (1.0 * n_deg) - for bl in degsh: - gf_to_symm[bl] << ss + for key in degsh: + if ss is None: + ss = gf_to_symm[key].copy() + ss.zero() + helper = ss.copy() + # get the transformation matrix + if isinstance(degsh, dict): + v, C = degsh[key] + else: + # for backward compatibility, allow degsh to be a list + v = numpy.eye(*ss.target_shape) + C = False + # the helper is in the basis where the blocks are all equal + helper.from_L_G_R(v.conjugate().transpose(), gf_to_symm[key], v) + if C: + helper << helper.transpose() + # average over all shells + ss += helper / (1.0 * n_deg) + # now put back the averaged gf to all shells + for key in degsh: + if isinstance(degsh, dict): + v, C = degsh[key] + else: + # for backward compatibility, allow degsh to be a list + v = numpy.eye(*ss.target_shape) + C = False + if C: + gf_to_symm[key].from_L_G_R(v, ss.transpose(), v.conjugate().transpose()) + else: + gf_to_symm[key].from_L_G_R(v, ss, v.conjugate().transpose()) def total_density(self, mu=None, iw_or_w="iw", with_Sigma=True, with_dc=True, broadening=None): r""" @@ -1610,3 +2053,9 @@ class SumkDFT(object): def __set_solver_to_sumk_block(self,value): self.block_structure.solver_to_sumk_block = value solver_to_sumk_block = property(__get_solver_to_sumk_block,__set_solver_to_sumk_block) + + def __get_deg_shells(self): + return self.block_structure.deg_shells + def __set_deg_shells(self,value): + self.block_structure.deg_shells = value + deg_shells = property(__get_deg_shells,__set_deg_shells) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 5ebe0867..b1dc6e8c 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -2,10 +2,10 @@ FILE(GLOB all_h5_files RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h5) file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/${all_h5_files} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) # Copy other files -FILE(COPY SrVO3.pmat SrVO3.struct SrVO3.outputs SrVO3.oubwin SrVO3.ctqmcout SrVO3.symqmc SrVO3.sympar SrVO3.parproj hk_convert_hamiltonian.hk LaVO3-Pnma_hr.dat LaVO3-Pnma.inp DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) +FILE(COPY SrVO3.pmat SrVO3.struct SrVO3.outputs SrVO3.oubwin SrVO3.ctqmcout SrVO3.symqmc SrVO3.sympar SrVO3.parproj SrIrO3_rot.h5 hk_convert_hamiltonian.hk LaVO3-Pnma_hr.dat LaVO3-Pnma.inp DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) # List all tests -set(all_tests wien2k_convert hk_convert w90_convert sumkdft_basic srvo3_Gloc srvo3_transp sigma_from_file blockstructure) +set(all_tests wien2k_convert hk_convert w90_convert sumkdft_basic srvo3_Gloc srvo3_transp sigma_from_file blockstructure analyze_block_structure_from_gf analyze_block_structure_from_gf2) foreach(t ${all_tests}) add_test(NAME ${t} COMMAND python ${CMAKE_CURRENT_SOURCE_DIR}/${t}.py) diff --git a/test/SrIrO3_rot.h5 b/test/SrIrO3_rot.h5 new file mode 100644 index 00000000..5ea92a78 Binary files /dev/null and b/test/SrIrO3_rot.h5 differ diff --git a/test/analyze_block_structure_from_gf.py b/test/analyze_block_structure_from_gf.py new file mode 100644 index 00000000..23ada79f --- /dev/null +++ b/test/analyze_block_structure_from_gf.py @@ -0,0 +1,232 @@ +from pytriqs.gf import * +from sumk_dft import SumkDFT +from scipy.linalg import expm +import numpy as np +from pytriqs.utility.comparison_tests import assert_gfs_are_close, assert_arrays_are_close, assert_block_gfs_are_close +from pytriqs.archive import * +import itertools + +# The full test checks all different possible combinations of conjugated +# blocks. This takes a few minutes. For a quick test, just checking one +# random value suffices. +# (this parameter affects the second test) +full_test = False + +####################################################################### +# First test # +# where we check the analyse_block_structure_from_gf function # +# for the SrIrO3_rot.h5 file # +####################################################################### + +beta = 40 +SK = SumkDFT(hdf_file = 'SrIrO3_rot.h5') +Sigma = SK.block_structure.create_gf(beta=beta) +SK.put_Sigma([Sigma]) +G = SK.extract_G_loc() + +# the original block structure +block_structure1 = SK.block_structure.copy() + +G_new = SK.analyse_block_structure_from_gf(G) + +# the new block structure +block_structure2 = SK.block_structure.copy() + +with HDFArchive('analyze_block_structure_from_gf.out.h5','w') as ar: + ar['bs1'] = block_structure1 + ar['bs2'] = block_structure2 + +# check whether the block structure is the same as in the reference +with HDFArchive('analyze_block_structure_from_gf.out.h5','r') as ar,\ + HDFArchive('analyze_block_structure_from_gf.ref.h5','r') as ar2: + assert ar['bs1'] == ar2['bs1'], 'bs1 not equal' + a1 = ar['bs2'] + a2 = ar2['bs2'] + assert a1==block_structure2, "writing/reading block structure incorrect" + # we set the deg_shells to None because the transformation matrices + # have a phase freedom and will, therefore, not be equal in general + a1.deg_shells = None + a2.deg_shells = None + assert a1==a2, 'bs2 not equal' + +# check if deg shells are correct +assert len(SK.deg_shells[0])==1, "wrong number of equivalent groups" + +# check if the Green's functions that are found to be equal in the +# routine are indeed equal +for d in SK.deg_shells[0]: + assert len(d)==2, "wrong number of shells in equivalent group" + # the convention is that for every degenerate shell, the transformation + # matrix v and the conjugate bool is saved + # then, + # maybe_conjugate1( v1^dagger G1 v1 ) = maybe_conjugate2( v2^dagger G2 v2 ) + # therefore, to test, we calculate + # maybe_conjugate( v^dagger G v ) + # for all degenerate shells and check that they are all equal + normalized_gfs = [] + for key in d: + normalized_gf = G_new[0][key].copy() + normalized_gf.from_L_G_R(d[key][0].conjugate().transpose(), G_new[0][key], d[key][0]) + if d[key][1]: + normalized_gf << normalized_gf.transpose() + normalized_gfs.append(normalized_gf) + for i in range(len(normalized_gfs)): + for j in range(i+1,len(normalized_gfs)): + assert_arrays_are_close(normalized_gfs[i].data, normalized_gfs[j].data, 1.e-5) + # the tails have to be compared using a relative error + for o in range(normalized_gfs[i].tail.order_min,normalized_gfs[i].tail.order_max+1): + if np.abs(normalized_gfs[i].tail[o][0,0]) < 1.e-10: + continue + assert np.max(np.abs((normalized_gfs[i].tail[o]-normalized_gfs[j].tail[o])/(normalized_gfs[i].tail[o][0,0]))) < 1.e-5, \ + "tails are different" + +####################################################################### +# Second test # +# where a Green's function is constructed from a random model # +# and the analyse_block_structure_from_gf function is tested for that # +# model # +####################################################################### + +# helper function to get random Hermitian matrix +def get_random_hermitian(dim): + herm = np.random.rand(dim,dim)+1.0j*np.random.rand(dim,dim) + herm = herm + herm.conjugate().transpose() + return herm + +# helper function to get random unitary matrix +def get_random_transformation(dim): + herm = get_random_hermitian(dim) + T = expm(1.0j*herm) + return T + +# we will conjugate the Green's function blocks according to the entries +# of conjugate_values +# for each of the 5 blocks that will be constructed, there is an entry +# True or False that says whether it will be conjugated +if full_test: + # in the full test we check all combinations + conjugate_values = list(itertools.product([False, True], repeat=5)) +else: + # in the quick test we check a random combination + conjugate_values = [np.random.rand(5)>0.5] + +for conjugate in conjugate_values: + # construct a random block-diagonal Hloc + Hloc = np.zeros((10,10), dtype=np.complex_) + # the Hloc of the first three 2x2 blocks is equal + Hloc0 = get_random_hermitian(2) + Hloc[:2,:2] = Hloc0 + Hloc[2:4,2:4] = Hloc0 + Hloc[4:6,4:6] = Hloc0 + # the Hloc of the last two 2x2 blocks is equal + Hloc1 = get_random_hermitian(2) + Hloc[6:8,6:8] = Hloc1 + Hloc[8:,8:] = Hloc1 + # construct the hybridization delta + # this is equal for all 2x2 blocks + V = get_random_hermitian(2) # the hopping elements from impurity to bath + b1 = np.random.rand() # the bath energy of the first bath level + b2 = np.random.rand() # the bath energy of the second bath level + delta = G[0]['ud'][:2,:2].copy() + delta[0,0] << (V[0,0]*V[0,0].conjugate()*inverse(Omega-b1)+V[0,1]*V[0,1].conjugate()*inverse(Omega-b2))/2.0 + delta[0,1] << (V[0,0]*V[1,0].conjugate()*inverse(Omega-b1)+V[0,1]*V[1,1].conjugate()*inverse(Omega-b2))/2.0 + delta[1,0] << (V[1,0]*V[0,0].conjugate()*inverse(Omega-b1)+V[1,1]*V[0,1].conjugate()*inverse(Omega-b2))/2.0 + delta[1,1] << (V[1,0]*V[1,0].conjugate()*inverse(Omega-b1)+V[1,1]*V[1,1].conjugate()*inverse(Omega-b2))/2.0 + # construct G + G[0].zero() + for i in range(0,10,2): + G[0]['ud'][i:i+2,i:i+2] << inverse(Omega-delta) + G[0]['ud'] << inverse(inverse(G[0]['ud']) - Hloc) + + # for testing symm_deg_gf below, we need this + # we construct it so that for every group of degenerate blocks of G[0], the + # mean of the blocks of G_noisy is equal to G[0] + G_noisy = G[0].copy() + noise1 = np.random.randn(*delta.target_shape) + G_noisy['ud'][:2,:2].data[:,:,:] += noise1 + G_noisy['ud'][2:4,2:4].data[:,:,:] -= noise1/2.0 + G_noisy['ud'][4:6,4:6].data[:,:,:] -= noise1/2.0 + noise2 = np.random.randn(*delta.target_shape) + G_noisy['ud'][6:8,6:8].data[:,:,:] += noise2 + G_noisy['ud'][8:,8:].data[:,:,:] -= noise2 + + # for testing backward-compatibility in symm_deg_gf, we need the + # un-transformed Green's functions + G_pre_transform = G[0].copy() + G_noisy_pre_transform = G_noisy.copy() + + # transform each block using a random transformation matrix + for i in range(0,10,2): + T = get_random_transformation(2) + G[0]['ud'][i:i+2,i:i+2].from_L_G_R(T, G[0]['ud'][i:i+2,i:i+2], T.conjugate().transpose()) + G_noisy['ud'][i:i+2,i:i+2].from_L_G_R(T, G_noisy['ud'][i:i+2,i:i+2], T.conjugate().transpose()) + # if that block shall be conjugated, go ahead and do it + if conjugate[i//2]: + G[0]['ud'][i:i+2,i:i+2] << G[0]['ud'][i:i+2,i:i+2].transpose() + G_noisy['ud'][i:i+2,i:i+2] << G_noisy['ud'][i:i+2,i:i+2].transpose() + + # analyse the block structure + G_new = SK.analyse_block_structure_from_gf(G, 1.e-7) + + # transform G_noisy etc. to the new block structure + G_noisy = SK.block_structure.convert_gf(G_noisy, block_structure1, beta = G_noisy.mesh.beta) + G_pre_transform = SK.block_structure.convert_gf(G_pre_transform, block_structure1, beta = G_noisy.mesh.beta) + G_noisy_pre_transform = SK.block_structure.convert_gf(G_noisy_pre_transform, block_structure1, beta = G_noisy.mesh.beta) + + assert len(SK.deg_shells[0]) == 2, "wrong number of equivalent groups found" + assert sorted([len(d) for d in SK.deg_shells[0]]) == [2,3], "wrong number of members in the equivalent groups found" + for d in SK.deg_shells[0]: + if len(d)==2: + assert 'ud_3' in d, "shell ud_3 missing" + assert 'ud_4' in d, "shell ud_4 missing" + if len(d)==3: + assert 'ud_0' in d, "shell ud_0 missing" + assert 'ud_1' in d, "shell ud_1 missing" + assert 'ud_2' in d, "shell ud_2 missing" + + # the convention is that for every degenerate shell, the transformation + # matrix v and the conjugate bool is saved + # then, + # maybe_conjugate1( v1^dagger G1 v1 ) = maybe_conjugate2( v2^dagger G2 v2 ) + # therefore, to test, we calculate + # maybe_conjugate( v^dagger G v ) + # for all degenerate shells and check that they are all equal + normalized_gfs = [] + for key in d: + normalized_gf = G_new[0][key].copy() + normalized_gf.from_L_G_R(d[key][0].conjugate().transpose(), G_new[0][key], d[key][0]) + if d[key][1]: + normalized_gf << normalized_gf.transpose() + normalized_gfs.append(normalized_gf) + for i in range(len(normalized_gfs)): + for j in range(i+1,len(normalized_gfs)): + # here, we use a threshold that is 1 order of magnitude less strict + # because of numerics + assert_gfs_are_close(normalized_gfs[i], normalized_gfs[j], 1.e-6) + + # now we check symm_deg_gf + # symmetrizing the GF has is has to leave it unchanged + G_new_symm = G_new[0].copy() + SK.symm_deg_gf(G_new_symm, 0) + assert_block_gfs_are_close(G_new[0], G_new_symm, 1.e-6) + + # symmetrizing the noisy GF, which was carefully constructed, + # has to give the same result as G_new[0] + SK.symm_deg_gf(G_noisy, 0) + assert_block_gfs_are_close(G_new[0], G_noisy, 1.e-6) + + # check backward compatibility of symm_deg_gf + # first, construct the old format of the deg shells + for ish in range(len(SK.deg_shells)): + for gr in range(len(SK.deg_shells[ish])): + SK.deg_shells[ish][gr] = SK.deg_shells[ish][gr].keys() + + # symmetrizing the GF as is has to leave it unchanged + G_new_symm << G_pre_transform + SK.symm_deg_gf(G_new_symm, 0) + assert_block_gfs_are_close(G_new_symm, G_pre_transform, 1.e-6) + + # symmetrizing the noisy GF pre transform, which was carefully constructed, + # has to give the same result as G_pre_transform + SK.symm_deg_gf(G_noisy_pre_transform, 0) + assert_block_gfs_are_close(G_noisy_pre_transform, G_pre_transform, 1.e-6) diff --git a/test/analyze_block_structure_from_gf.ref.h5 b/test/analyze_block_structure_from_gf.ref.h5 new file mode 100644 index 00000000..7acea34d Binary files /dev/null and b/test/analyze_block_structure_from_gf.ref.h5 differ diff --git a/test/analyze_block_structure_from_gf2.py b/test/analyze_block_structure_from_gf2.py new file mode 100644 index 00000000..d371a9c5 --- /dev/null +++ b/test/analyze_block_structure_from_gf2.py @@ -0,0 +1,115 @@ +from pytriqs.gf import * +from sumk_dft import SumkDFT +import numpy as np +from pytriqs.utility.comparison_tests import assert_block_gfs_are_close + +# here we test the SK.analyze_block_structure_from_gf function +# with GfReFreq, GfReTime + + +# helper function to get random Hermitian matrix +def get_random_hermitian(dim): + herm = np.random.rand(dim,dim)+1.0j*np.random.rand(dim,dim) + herm = herm + herm.conjugate().transpose() + return herm + +# helper function to get random unitary matrix +def get_random_transformation(dim): + herm = get_random_hermitian(dim) + T = expm(1.0j*herm) + return T + +# construct a random block-diagonal Hloc +Hloc = np.zeros((10,10), dtype=np.complex_) +# the Hloc of the first three 2x2 blocks is equal +Hloc0 = get_random_hermitian(2) +Hloc[:2,:2] = Hloc0 +Hloc[2:4,2:4] = Hloc0 +Hloc[4:6,4:6] = Hloc0 +# the Hloc of the last two 2x2 blocks is equal +Hloc1 = get_random_hermitian(2) +Hloc[6:8,6:8] = Hloc1 +Hloc[8:,8:] = Hloc1 +# construct the hybridization delta +# this is equal for all 2x2 blocks +V = get_random_hermitian(2) # the hopping elements from impurity to bath +b1 = np.random.rand() # the bath energy of the first bath level +b2 = np.random.rand() # the bath energy of the second bath level +delta = GfReFreq(window=(-5,5), indices=range(2), n_points=1001) +delta[0,0] << (V[0,0]*V[0,0].conjugate()*inverse(Omega-b1)+V[0,1]*V[0,1].conjugate()*inverse(Omega-b2+0.02j))/2.0 +delta[0,1] << (V[0,0]*V[1,0].conjugate()*inverse(Omega-b1)+V[0,1]*V[1,1].conjugate()*inverse(Omega-b2+0.02j))/2.0 +delta[1,0] << (V[1,0]*V[0,0].conjugate()*inverse(Omega-b1)+V[1,1]*V[0,1].conjugate()*inverse(Omega-b2+0.02j))/2.0 +delta[1,1] << (V[1,0]*V[1,0].conjugate()*inverse(Omega-b1)+V[1,1]*V[1,1].conjugate()*inverse(Omega-b2+0.02j))/2.0 +# construct G +G = BlockGf(name_block_generator=(('ud',GfReFreq(window=(-5,5), indices=range(10), n_points=1001)),)) +for i in range(0,10,2): + G['ud'][i:i+2,i:i+2] << inverse(Omega-delta+0.02j) +G['ud'] << inverse(inverse(G['ud']) - Hloc) + + +SK = SumkDFT(hdf_file = 'SrIrO3_rot.h5', use_dft_blocks=False) +G_new = SK.analyse_block_structure_from_gf([G]) +G_new_symm = G_new[0].copy() +SK.symm_deg_gf(G_new_symm, 0) +assert_block_gfs_are_close(G_new[0], G_new_symm) + + +assert SK.gf_struct_sumk == [[('ud', [0, 1, 2, 3, 4, 5, 6, 7, 8, 9])], [('ud', [0, 1, 2, 3, 4, 5, 6, 7, 8, 9])]],\ + "wrong gf_struct_sumk" +for i in range(5): + assert 'ud_{}'.format(i) in SK.gf_struct_solver[0], "missing block" + assert SK.gf_struct_solver[0]['ud_{}'.format(i)] == range(2), "wrong block size" +for i in range(10): + assert SK.sumk_to_solver[0]['ud',i] == ('ud_{}'.format(i/2), i%2), "wrong mapping" + +assert len(SK.deg_shells[0]) == 2, "wrong number of equivalent groups found" +assert sorted([len(d) for d in SK.deg_shells[0]]) == [2,3], "wrong number of members in the equivalent groups found" +for d in SK.deg_shells[0]: + if len(d)==2: + assert 'ud_3' in d, "shell ud_3 missing" + assert 'ud_4' in d, "shell ud_4 missing" + if len(d)==3: + assert 'ud_0' in d, "shell ud_0 missing" + assert 'ud_1' in d, "shell ud_1 missing" + assert 'ud_2' in d, "shell ud_2 missing" + + + +def get_delta_from_mesh(mesh): + w0 = None + for w in mesh: + if w0 is None: + w0 = w + else: + return w-w0 + +Gt = BlockGf(name_block_generator = [(name, + GfReTime(window=(-np.pi*(len(block.mesh)-1) / (len(block.mesh)*get_delta_from_mesh(block.mesh)), np.pi*(len(block.mesh)-1) / (len(block.mesh)*get_delta_from_mesh(block.mesh))), + n_points=len(block.mesh), + indices=block.indices)) for name, block in G]) + +Gt['ud'].set_from_inverse_fourier(G['ud']) + +G_new = SK.analyse_block_structure_from_gf([Gt]) +G_new_symm = G_new[0].copy() +SK.symm_deg_gf(G_new_symm, 0) +assert_block_gfs_are_close(G_new[0], G_new_symm) + +assert SK.gf_struct_sumk == [[('ud', [0, 1, 2, 3, 4, 5, 6, 7, 8, 9])], [('ud', [0, 1, 2, 3, 4, 5, 6, 7, 8, 9])]],\ + "wrong gf_struct_sumk" +for i in range(5): + assert 'ud_{}'.format(i) in SK.gf_struct_solver[0], "missing block" + assert SK.gf_struct_solver[0]['ud_{}'.format(i)] == range(2), "wrong block size" +for i in range(10): + assert SK.sumk_to_solver[0]['ud',i] == ('ud_{}'.format(i/2), i%2), "wrong mapping" + +assert len(SK.deg_shells[0]) == 2, "wrong number of equivalent groups found" +assert sorted([len(d) for d in SK.deg_shells[0]]) == [2,3], "wrong number of members in the equivalent groups found" +for d in SK.deg_shells[0]: + if len(d)==2: + assert 'ud_3' in d, "shell ud_3 missing" + assert 'ud_4' in d, "shell ud_4 missing" + if len(d)==3: + assert 'ud_0' in d, "shell ud_0 missing" + assert 'ud_1' in d, "shell ud_1 missing" + assert 'ud_2' in d, "shell ud_2 missing" diff --git a/test/blockstructure.in.h5 b/test/blockstructure.in.h5 index dccb8ef0..6e1bb469 100644 Binary files a/test/blockstructure.in.h5 and b/test/blockstructure.in.h5 differ diff --git a/test/blockstructure.py b/test/blockstructure.py index 33f1a570..ea2c7310 100644 --- a/test/blockstructure.py +++ b/test/blockstructure.py @@ -21,7 +21,8 @@ sk_pick1 = BlockStructure(gf_struct_sumk = SK.gf_struct_sumk, gf_struct_solver = SK.gf_struct_solver, solver_to_sumk = SK.solver_to_sumk, sumk_to_solver = SK.sumk_to_solver, - solver_to_sumk_block = SK.solver_to_sumk_block) + solver_to_sumk_block = SK.solver_to_sumk_block, + deg_shells = SK.deg_shells) assert sk_pick1 == pick1, 'constructing block structure from SumkDFT properties failed' # check pick_gf_struct_sumk diff --git a/test/blockstructure.ref.h5 b/test/blockstructure.ref.h5 index b290411c..c9eb4230 100644 Binary files a/test/blockstructure.ref.h5 and b/test/blockstructure.ref.h5 differ