Merge remote-tracking branch 'gernot/analyze_block_structure_from_gf' into analyze_block_structure_from_gf

python/block_structure.py

@@ -2,6 +2,7 @@ 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):
@@ -40,12 +41,14 @@ class BlockStructure(object):
                       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,7 +102,8 @@ 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.
@@ -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
+        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,14 +324,16 @@ 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]
@@ -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
@@ -376,7 +392,8 @@ class BlockStructure(object):
             return False
 
         for prop in [ "gf_struct_sumk", "gf_struct_solver",
-                "solver_to_sumk", "sumk_to_solver", "solver_to_sumk_block"]:
+                "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
@@ -389,7 +406,7 @@ class BlockStructure(object):
 
         ret = {}
         for element in [ "gf_struct_sumk", "gf_struct_solver",
-                         "solver_to_sumk_block"]:
+                         "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

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 <dft.block_structure.BlockStructure>`
+        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 <dft.block_structure.BlockStructure>`
+        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)

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)

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)

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"

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