BlockStructure class for manipulating GF structures

BlockStructure stores gf_struct_solver and gf_struct_sumk
and the mappings between them. It allows to modify it, and
save it to HDF (see issue #67).

doc/documentation.rst

@@ -45,6 +45,7 @@ This is the reference manual for the python routines.
    reference/sumk_dft_tools
    reference/symmetry
    reference/transbasis
+   reference/block_structure
 
 
 FAQs

doc/reference/block_structure.rst

@@ -0,0 +1,22 @@
+Block Structure
+===============
+
+The `BlockStructure` class allows to change and manipulate
+Green's functions structures and mappings from sumk to solver.
+
+The block structure can also be written to and read from HDF files.
+
+.. warning::
+
+    Do not write the individual elements of this class to a HDF file,
+    as they belong together and changing one without the other can 
+    result in unexpected results. Always write the BlockStructure 
+    object as a whole.
+
+    Writing the sumk_to_solver and solver_to_sumk elements
+    individually is not implemented.
+
+.. autoclass:: dft.block_structure.BlockStructure
+   :members:
+   :show-inheritance:
+   

python/__init__.py

@@ -22,8 +22,9 @@
 
 from sumk_dft import SumkDFT
 from symmetry import Symmetry
+from block_structure import BlockStructure
 from sumk_dft_tools import SumkDFTTools
 from converters import *
 
 __all__ = ['SumkDFT', 'Symmetry', 'SumkDFTTools',
-           'Wien2kConverter', 'HkConverter']
+           'Wien2kConverter', 'HkConverter','BlockStructure']

python/block_structure.py

@@ -0,0 +1,442 @@
+import copy
+import numpy as np
+from pytriqs.gf.local import GfImFreq, BlockGf
+from ast import literal_eval
+from warnings import warn
+
+class BlockStructure(object):
+    """ Contains information about the Green function structure.
+
+    This class contains information about the structure of the solver 
+    and sumk Green functions and the mapping between them.
+
+    Parameters
+    ----------
+    gf_struct_sumk : list of list of tuple
+        gf_struct_sumk[ish][idx] = (block_name,list of indices in block)
+
+        for correlated shell ish; idx is just a counter in the list
+    gf_struct_solver : list of dict
+        gf_struct_solver[ish][block] = list of indices in that block
+
+        for *inequivalent* correlated shell ish
+    solver_to_sumk : list of dict
+        solver_to_sumk[ish][(from_block,from_idx)] = (to_block,to_idx)
+
+        maps from the solver block and index to the sumk block and index
+        for *inequivalent* correlated shell ish
+    sumk_to_solver : list of dict
+        sumk_to_solver[ish][(from_block,from_idx)] = (to_block,to_idx)
+
+        maps from the sumk block and index to the solver block and index
+        for *inequivalent* correlated shell ish
+    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 
+        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):
+        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
+
+    @classmethod
+    def full_structure(cls,gf_struct,corr_to_inequiv):
+        """ Construct structure that maps to itself.
+
+        This has the same structure for sumk and solver, and the
+        mapping solver_to_sumk and sumk_to_solver is one-to-one.
+
+        Parameters
+        ----------
+        gf_struct : list of dict
+            gf_struct[ish][block] = list of indices in that block
+
+            for (inequivalent) correlated shell ish
+        corr_to_inequiv : list
+            gives the mapping from correlated shell csh to inequivalent
+            correlated shell icsh, so that corr_to_inequiv[csh]=icsh
+            e.g. SumkDFT.corr_to_inequiv
+
+            if None, each inequivalent correlated shell is supposed to
+            be correspond to just one correlated shell with the same
+            index; there is not default, None has to be set explicitly!
+        """
+
+        solver_to_sumk = []
+        s2sblock = []
+        gs_sumk = []
+        for ish in range(len(gf_struct)):
+            so2su = {}
+            so2sublock = {}
+            gss = []
+            for block in gf_struct[ish]:
+                so2sublock[block]=block
+                for ind in gf_struct[ish][block]:
+                    so2su[(block,ind)]=(block,ind)
+                gss.append((block,gf_struct[ish][block]))
+            solver_to_sumk.append(so2su)
+            s2sblock.append(so2sublock)
+            gs_sumk.append(gss)
+
+        # gf_struct_sumk is not given for each inequivalent correlated
+        # shell, but for every correlated shell!
+        if corr_to_inequiv is not None:
+            gs_sumk_all = [None]*len(corr_to_inequiv)
+            for i in range(len(corr_to_inequiv)):
+                gs_sumk_all[i] = gs_sumk[corr_to_inequiv[i]]
+        else:
+            gs_sumk_all = gs_sumk
+
+        return cls(gf_struct_solver=copy.deepcopy(gf_struct),
+                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)
+
+    def pick_gf_struct_solver(self,new_gf_struct):
+        """ 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 
+        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: 
+        
+        [{'up':[0,1],'down':[0,1],'left':[0,1]}]
+
+        to choose the 0th index of the up block and the 1st index of
+        the down block and drop the left block, the new_gf_struct would
+        have to be
+
+        [{'up':[0],'down':[1]}]
+
+        Note that the indices will be renamed to be a 0-based
+        sequence of integers, i.e. the new structure will actually
+        be  [{'up':[0],'down':[0]}].
+
+        For dropped indices, sumk_to_solver will map to (None,None).
+
+        Parameters
+        ----------
+        new_gf_struct : list of dict
+            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]
+
+            # create new solver_to_sumk
+            so2su={}
+            so2su_block = {}
+            for blk,idxs in gf_struct.items():
+                for i in range(len(idxs)):
+                    so2su[(blk,i)]=self.solver_to_sumk[ish][(blk,idxs[i])]
+                    so2su_block[blk]=so2su[(blk,i)][0]
+            self.solver_to_sumk[ish] = so2su
+            self.solver_to_sumk_block[ish] = so2su_block
+            # create new sumk_to_solver
+            for k,v in self.sumk_to_solver[ish].items():
+                blk,ind=v
+                if blk in gf_struct and ind in gf_struct[blk]:
+                    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) 
+            # 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. 
+
+        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 
+        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: 
+        
+        [{'up':[0,1],'down':[0,1],'left':[0,1]}]
+
+        to choose the 0th index of the up block and the 1st index of
+        the down block and drop the left block, the new_gf_struct would
+        have to be
+
+        [{'up':[0],'down':[1]}]
+
+        Note that the indices will be renamed to be a 0-based
+        sequence of integers.
+
+        For dropped indices, sumk_to_solver will map to (None,None).
+
+        Parameters
+        ----------
+        new_gf_struct : list of dict
+            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 
+            but the sumk Gf.
+        """
+
+
+        gfs = []
+        # construct gfs, which is the equivalent of new_gf_struct
+        # but according to the solver Gf, by using the sumk_to_solver
+        # mapping
+        for ish in range(len(new_gf_struct)):
+            gfs.append({})
+            for block in new_gf_struct[ish].keys():
+                for ind in new_gf_struct[ish][block]:
+                    ind_sol = self.sumk_to_solver[ish][(block,ind)]
+                    if not ind_sol[0] in gfs[ish]:
+                        gfs[ish][ind_sol[0]]=[]
+                    gfs[ish][ind_sol[0]].append(ind_sol[1])
+        self.pick_gf_struct_solver(gfs)
+
+
+    def map_gf_struct_solver(self,mapping):
+        """ Map the Green function structure from one struct to another.
+
+        Parameters
+        ----------
+        mapping : list of dict
+            the dict consists of elements 
+            (from_block,from_index) : (to_block,to_index)
+            that maps from one structure to the other
+        """
+
+        for ish in range(len(mapping)):
+            gf_struct = {}
+            so2su = {}
+            su2so = {}
+            so2su_block = {}
+            for frm,to in mapping[ish].iteritems():
+                if not to[0] in gf_struct:
+                    gf_struct[to[0]]=[]
+                gf_struct[to[0]].append(to[1])
+
+                so2su[to]=self.solver_to_sumk[ish][frm]
+                su2so[self.solver_to_sumk[ish][frm]]=to
+                if to[0] in so2su_block:
+                    if so2su_block[to[0]] != \
+                        self.solver_to_sumk_block[ish][frm[0]]:
+                            warn("solver block '{}' maps to more than one sumk block: '{}', '{}'".format(
+                                    to[0],so2su_block[to[0]],self.solver_to_sumk_block[ish][frm[0]]))
+                else:
+                    so2su_block[to[0]]=\
+                        self.solver_to_sumk_block[ish][frm[0]]
+            for k in self.sumk_to_solver[ish].keys():
+                if not k in su2so:
+                    su2so[k] = (None,None)
+            self.gf_struct_solver[ish]=gf_struct
+            self.solver_to_sumk[ish]=so2su
+            self.sumk_to_solver[ish]=su2so
+            self.solver_to_sumk_block[ish]=so2su_block
+
+    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 
+        supply beta as keyword argument.
+
+        Parameters
+        ----------
+        ish : int
+            shell index
+        gf_function : constructor
+            function used to construct the Gf objects constituting the
+            individual blocks; default: GfImFreq
+        **kwargs :
+            options passed on to the Gf constructor for the individual
+            blocks
+        """
+
+        names = self.gf_struct_solver[ish].keys()
+        blocks=[]
+        for n in names:
+            G = gf_function(indices=self.gf_struct_solver[ish][n],**kwargs)
+            blocks.append(G)
+        G = BlockGf(name_list = names, block_list = blocks)
+        return G
+
+
+    def convert_gf(self,G,G_struct,ish=0,show_warnings=True,**kwargs):
+        """ Convert BlockGf from its structure to this structure.
+
+        .. warning::
+
+            Elements that are zero in the new structure due to
+            the new block structure will be just ignored, thus 
+            approximated to zero.
+
+        Parameters
+        ----------
+        G : BlockGf
+            the Gf that should be converted
+        G_struct : GfStructure
+            the structure ofthat G 
+        ish : int
+            shell index
+        show_warnings : bool
+            whether to show warnings when elements of the Green's 
+            function get thrown away
+        **kwargs :
+            options passed to the constructor for the new Gf
+        """
+        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]:
+                for i2 in G_struct.gf_struct_solver[ish][block]:
+                    i1_sumk = G_struct.solver_to_sumk[ish][(block,i1)]
+                    i2_sumk = G_struct.solver_to_sumk[ish][(block,i2)]
+                    i1_sol = self.sumk_to_solver[ish][i1_sumk]
+                    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))
+                        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))
+                        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. 
+        
+        .. warning::
+
+            In general, this will throw away non-zero elements of the
+            Green's function. Be sure to verify whether this approximation
+            is justified.
+        """
+
+        self.gf_struct_solver=[]
+        self.solver_to_sumk=[]
+        self.solver_to_sumk_block=[]
+        for ish in range(len(self.sumk_to_solver)):
+            self.gf_struct_solver.append({})
+            self.solver_to_sumk.append({})
+            self.solver_to_sumk_block.append({})
+            for frm,to in self.sumk_to_solver[ish].iteritems():
+                if to[0] is not None:
+                    self.gf_struct_solver[ish][frm[0]+'_'+str(frm[1])]=[0]
+                    self.sumk_to_solver[ish][frm]=(frm[0]+'_'+str(frm[1]),0)
+                    self.solver_to_sumk[ish][(frm[0]+'_'+str(frm[1]),0)]=frm
+                    self.solver_to_sumk_block[ish][frm[0]+'_'+str(frm[1])]=frm[0]
+
+    def __eq__(self,other):
+        def compare(one,two):
+            if type(one)!=type(two):
+                return False
+            if one is None and two is None:
+                return True
+            if isinstance(one,list) or isinstance(one,tuple):
+                if len(one) != len(two):
+                    return False
+                for x,y in zip(one,two):
+                    if not compare(x,y):
+                        return False
+                return True
+            elif isinstance(one,int): 
+                return one==two
+            elif isinstance(one,str):
+                return one==two
+            elif isinstance(one,dict):
+                if set(one.keys()) != set(two.keys()):
+                    return False
+                for k in set(one.keys()).intersection(two.keys()):
+                    if not compare(one[k],two[k]):
+                        return False
+                return True
+            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"]:
+            if not compare(getattr(self,prop),getattr(other,prop)):
+                return False
+        return True
+
+    def copy(self):
+        return copy.deepcopy(self)
+
+    def __reduce_to_dict__(self):
+        """ Reduce to dict for HDF5 export."""
+
+        ret = {}
+        for element in [ "gf_struct_sumk", "gf_struct_solver", 
+                         "solver_to_sumk_block"]:
+            ret[element] = getattr(self,element)
+
+        def construct_mapping(mapping):
+            d = []
+            for ish in range(len(mapping)):
+                d.append({})
+                for k,v in mapping[ish].iteritems():
+                    d[ish][repr(k)] = repr(v)
+            return d
+
+        ret['solver_to_sumk']=construct_mapping(self.solver_to_sumk)
+        ret['sumk_to_solver']=construct_mapping(self.sumk_to_solver)
+        return ret
+
+    @classmethod
+    def __factory_from_dict__(cls,name,D) :
+        """ Create from dict for HDF5 import."""
+
+        def reconstruct_mapping(mapping):
+            d = []
+            for ish in range(len(mapping)):
+                d.append({})
+                for k,v in mapping[ish].iteritems():
+                    # literal_eval is a saje alternative to eval
+                    d[ish][literal_eval(k)] = literal_eval(v)
+            return d
+
+        D['solver_to_sumk']=reconstruct_mapping(D['solver_to_sumk'])
+        D['sumk_to_solver']=reconstruct_mapping(D['sumk_to_solver'])
+        return cls(**D)
+
+    def __str__(self):
+        s=''
+        s+= "gf_struct_sumk "+str( self.gf_struct_sumk)+'\n'
+        s+= "gf_struct_solver "+str(self.gf_struct_solver)+'\n'
+        s+= "solver_to_sumk_block "+str(self.solver_to_sumk_block)+'\n'
+        for el in ['solver_to_sumk','sumk_to_solver']:
+            s+=el+'\n'
+            element=getattr(self,el)
+            for ish in range(len(element)):
+                s+=' shell '+str(ish)+'\n'
+                def keyfun(el):
+                    return '{}_{:05d}'.format(el[0],el[1])
+                keys = sorted(element[ish].keys(),key=keyfun)
+                for k in keys:
+                    s+='  '+str(k)+str(element[ish][k])+'\n'
+        return s
+
+from pytriqs.archive.hdf_archive_schemes import register_class
+register_class(BlockStructure)

python/sumk_dft.py

@@ -27,12 +27,13 @@ from pytriqs.gf.local import *
 import pytriqs.utility.mpi as mpi
 from pytriqs.archive import *
 from symmetry import *
+from block_structure import BlockStructure
 from sets import Set
 from itertools import product
 from warnings import warn
 
 
-class SumkDFT:
+class SumkDFT(object):
     """This class provides a general SumK method for combining ab-initio code and pytriqs."""
 
     def __init__(self, hdf_file, h_field=0.0, use_dft_blocks=False,
@@ -53,6 +54,9 @@ class SumkDFT:
         use_dft_blocks : boolean, optional
                          If True, the local Green's function matrix for each spin is divided into smaller blocks 
                           with the block structure determined from the DFT density matrix of the corresponding correlated shell.
+
+                         Alternatively and additionally, the block structure can be analyzed using :meth:`analyse_block_structure <dft.sumk_dft.SumkDFT.analyse_block_structure>`
+                         and manipulated using the SumkDFT.block_structre attribute (see :class:`BlockStructure <dft.block_structure.BlockStructure>`).
         dft_data : string, optional
                    Name of hdf5 subgroup in which DFT data for projector and lattice Green's function construction are stored.
         symmcorr_data : string, optional
@@ -112,6 +116,8 @@ class SumkDFT:
                     self.spin_names_to_ind[iso][
                         self.spin_block_names[iso][isp]] = isp * self.SP
 
+            self.block_structure = BlockStructure()
+
             # GF structure used for the local things in the k sums
             # Most general form allowing for all hybridisation, i.e. largest
             # blocks possible
@@ -221,6 +227,9 @@ class SumkDFT:
         if not subgrp in ar:
             ar.create_group(subgrp)
         for it in things_to_save:
+            if it in [ "gf_struct_sumk", "gf_struct_solver",
+                    "solver_to_sumk", "sumk_to_solver", "solver_to_sumk_block"]:
+                warn("It is not recommended to save '{}' individually. Save 'block_structure' instead.".format(it))
             try:
                 ar[subgrp][it] = getattr(self, it)
             except:
@@ -719,7 +728,7 @@ class SumkDFT:
         r"""
         Determines the block structure of local Green's functions by analysing the structure of 
         the corresponding density matrices and the local Hamiltonian. The resulting block structures 
-        for correlated shells are stored in self.gf_struct_solver list.
+        for correlated shells are stored in the :class:`SumkDFT.block_structure <dft.block_structure.BlockStructure>` attribute.
 
         Parameters
         ----------
@@ -1505,3 +1514,35 @@ class SumkDFT:
         atomlst = [shells[i]['atom'] for i in range(len(shells))]
         n_atoms = len(set(atomlst))
         return n_atoms
+
+    # The following methods are here to ensure backward-compatibility
+    # after introducing the block_structure class
+    def __get_gf_struct_sumk(self):
+        return self.block_structure.gf_struct_sumk
+    def __set_gf_struct_sumk(self,value):
+        self.block_structure.gf_struct_sumk = value
+    gf_struct_sumk = property(__get_gf_struct_sumk,__set_gf_struct_sumk)
+
+    def __get_gf_struct_solver(self):
+        return self.block_structure.gf_struct_solver
+    def __set_gf_struct_solver(self,value):
+        self.block_structure.gf_struct_solver = value
+    gf_struct_solver = property(__get_gf_struct_solver,__set_gf_struct_solver)
+
+    def __get_solver_to_sumk(self):
+        return self.block_structure.solver_to_sumk
+    def __set_solver_to_sumk(self,value):
+        self.block_structure.solver_to_sumk = value
+    solver_to_sumk = property(__get_solver_to_sumk,__set_solver_to_sumk)
+
+    def __get_sumk_to_solver(self):
+        return self.block_structure.sumk_to_solver
+    def __set_sumk_to_solver(self,value):
+        self.block_structure.sumk_to_solver = value
+    sumk_to_solver = property(__get_sumk_to_solver,__set_sumk_to_solver)
+
+    def __get_solver_to_sumk_block(self):
+        return self.block_structure.solver_to_sumk_block
+    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)

test/CMakeLists.txt

@@ -14,3 +14,4 @@ triqs_add_python_test(sumkdft_basic)
 triqs_add_python_test(srvo3_Gloc)
 triqs_add_python_test(srvo3_transp)
 triqs_add_python_test(sigma_from_file)
+triqs_add_python_test(blockstructure)

test/blockstructure.py

@@ -0,0 +1,83 @@
+from pytriqs.applications.dft.sumk_dft import *
+from pytriqs.utility.h5diff import h5diff
+from pytriqs.gf.local import *
+from pytriqs.utility.comparison_tests import assert_block_gfs_are_close
+from pytriqs.applications.dft import BlockStructure
+
+SK = SumkDFT('blockstructure.in.h5',use_dft_blocks=True)
+
+original_bs = SK.block_structure
+
+# check pick_gf_struct_solver
+pick1 = original_bs.copy()
+pick1.pick_gf_struct_solver([{'up_0': [1], 'up_1': [0], 'down_1': [0]}])
+
+# check loading a block_structure from file
+SK.block_structure = SK.load(['block_structure'],'mod')[0]
+assert SK.block_structure == pick1, 'loading SK block structure from file failed'
+
+# check SumkDFT backward compatibility
+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)
+assert sk_pick1 == pick1, 'constructing block structure from SumkDFT properties failed'
+
+# check pick_gf_struct_sumk
+pick2 = original_bs.copy()
+pick2.pick_gf_struct_sumk([{'up': [1, 2], 'down': [0,1]}])
+
+# check map_gf_struct_solver
+mapping = [{ ('down_0', 0):('down', 0),
+             ('down_0', 1):('down', 2),
+             ('down_1', 0):('down', 1),
+             ('up_0', 0)  :('down_1', 0),
+             ('up_0', 1)  :('up_0', 0) }]
+map1 = original_bs.copy()
+map1.map_gf_struct_solver(mapping)
+
+# check create_gf
+G1 = original_bs.create_gf(beta=40,n_points=3)
+i = 1
+for block,gf in G1:
+    gf << SemiCircular(i)
+    i+=1
+
+# check approximate_as_diagonal
+offd = original_bs.copy()
+offd.approximate_as_diagonal()
+
+# check map_gf_struct_solver
+G2 = map1.convert_gf(G1,original_bs,beta=40,n_points=3,show_warnings=False)
+
+# check full_structure
+full = BlockStructure.full_structure([{'up_0': [0, 1], 'up_1': [0], 'down_1': [0], 'down_0': [0, 1]}],None)
+
+# check __eq__
+assert full==full, 'equality not correct (equal structures not equal)'
+assert pick1==pick1, 'equality not correct (equal structures not equal)'
+assert pick1!=pick2, 'equality not correct (different structures not different)'
+assert original_bs!=offd, 'equality not correct (different structures not different)'
+
+if mpi.is_master_node():
+    with HDFArchive('blockstructure.out.h5','w') as ar:
+        ar['original_bs'] = original_bs
+        ar['pick1'] = pick1
+        ar['pick2'] = pick2
+        ar['map1'] = map1
+        ar['offd'] = offd
+        ar['G1'] = G1
+        ar['G2'] = G2
+        ar['full'] = full
+
+    # cannot use h5diff because BlockStructure testing is not implemented
+    # there (and seems difficult to implement because it would mix triqs
+    # and dft_tools)
+    with HDFArchive('blockstructure.out.h5','r') as ar,\
+         HDFArchive('blockstructure.ref.h5','r') as ar2:
+            for k in ar2:
+                if isinstance(ar[k],BlockGf):
+                    assert_block_gfs_are_close(ar[k],ar2[k],1.e-6)
+                else:
+                    assert ar[k]==ar2[k], '{} not equal'.format(k)