Merge branch 'gernotsandmychanges' into HEAD

doc/documentation.rst

@@ -33,6 +33,8 @@ DFT+DMFT
 
    guide/dftdmft_singleshot
    guide/dftdmft_selfcons
+   guide/Sr2RuO4
+   guide/BasisRotation
 
 Postprocessing
 --------------

doc/guide/BasisRotation.rst

@@ -0,0 +1,76 @@
+.. _plovasp:
+
+Numerical Treatment of the Sign-Problem: Basis Rotations
+=======
+
+When performing calculations with off-diagonal complex hybridisation or local Hamiltonian, one is
+often limited by the fermionic sign-problem. However, as the sign is no
+physical observable, one can try to improve the calculation by rotating
+to another basis.
+
+While the choice of basis to perform the calculation in can be chosen
+arbitrarily, two choices which have shown good results are the basis
+which diagonalizes the local Hamiltonian or the density matrix of then
+system.
+
+The transformation matrix can be stored in the :class:`BlockStructure` and the
+transformation is automatically performed when using :class:`SumkDFT`'s :meth:`extract_G_loc`
+and :meth:`put_Sigma` (see below).
+
+
+Finding the Transformation Matrix
+-----------------
+
+The :class:`TransBasis` class offers a simple mehod to calculate the transformation
+matrices to a basis where either the local Hamiltonian or the density matrix
+is diagonal::
+
+    from triqs_dft_tools.trans_basis import TransBasis
+    TB = TransBasis(SK)
+    TB.calculate_diagonalisation_matrix(prop_to_be_diagonal='eal', calc_in_solver_blocks = True)
+
+    SK.block_structure.transformation = [{'ud':TB.w}]
+
+
+
+Transforming Green's functions manually
+-----------------
+
+One can transform Green's functions manually using the :meth:`convert_gf` method::
+
+    # Rotate a Green's function from solver-space to sumk-space
+    new_gf = block_structure.convert_gf(old_gf, space_from='solver', space_to='sumk')
+
+
+
+Automatic transformation during the DMFT loop
+-----------------
+
+During a DMFT loop one is switching back and forth between Sumk-Space and Solver-Space
+in each iteration. Once the block_structure.transformation property is set, this can be
+done automatically::
+
+    for it in range(iteration_offset, iteration_offset + n_iterations):
+        # every GF is in solver space here
+        S.G0_iw << inverse(S.Sigma_iw + inverse(S.G_iw))
+
+        # solve the impurity in solver space -> hopefully better sign
+        S.solve(h_int = H, **p)
+
+        # calc_dc(..., transform = True) by default
+        SK.calc_dc(S.G_iw.density(), U_interact=U, J_hund=J, orb=0, use_dc_formula=DC_type)
+        
+        # put_Sigma(..., transform_to_sumk_blocks = True) by default
+        SK.put_Sigma([S.Sigma_iw])
+        
+        SK.calc_mu()
+
+        # extract_G_loc(..., transform_to_solver_blocks = True) by default
+        S.G_iw << SK.extract_G_loc()[0]
+
+.. warning::
+  One must not forget to also transform the interaction Hamiltonian to the diagonal basis!
+  This can be done with the :meth:`transform_U_matrix` method. However, due to different 
+  conventions in this method, one must pass the conjugated version of the transformation matrix::
+  
+  U_trans = transform_U_matrix(U, SK.block_structure.transformation[0]['ud'].conjugate())

doc/guide/Sr2RuO4.rst

@@ -0,0 +1,40 @@
+.. _Sr2RuO4:
+
+Spin-orbit coupled calculations (single-shot)
+=============================================
+
+There are two main ways of including the spin-orbit coupling (SO) term into
+DFT+DMFT calculations:
+
+-   by performing a DFT calculation including SO and then doing a DMFT calculation on top, or
+-   by performing a DFT calculation without SO and then adding the SO term on the model level.
+
+Treatment of SO in DFT
+----------------------
+
+For now, TRIQS/DFTTools does only work with Wien2k when performing calculations with SO.
+Of course, the general Hk framework is also possible.
+But the way VASP treats SO is fundamentally different to the way Wien2k treats it and the interface does not handle that at the moment.
+
+Therefore, this guide assumes that Wien2k is being used.
+
+First, a Wien2k calculation including SO has to be performed.
+For details, we refer the reader to the documentation of Wien2k.
+The interface to Wien2k only works when the DFT calculation is done both spin-polarized and with SO (that means that you have to initialize the Wien2k calculation accordingly and then run with ``runsp -sp``).
+
+Performing the projection
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Note that the final ``x lapw2 -almd -so -up`` and ``x lapw2 -almd -so -dn`` have to be run *on a single core*, which implies that, before, ``x lapw1 -up``, ``x lapw2 -dn``, and ``x lapwso -up`` have to be run in single-core mode (once).
+
+In the ``case.indmftpr`` file, the spin-orbit flag has to be set to ``1`` for the correlated atoms.
+For example, for the compound Sr\ :sub:`2`\ RuO\ :sub:`4`, with the struct file :download:`Sr2RuO4.struct <Sr2RuO4/Sr2RuO4.struct>`, we would e.g. use the ``indmftpr`` file :download:`found here <Sr2RuO4/Sr2RuO4.indmftpr>`.
+Then, ``dmftproj -sp -so`` has to be called.
+As usual, it is important to check for warnings (e.g., about eigenvalues of the overlap matrix) in the output of ``dmftproj`` and adapt the window until these warnings disappear.
+
+Note that in presence of SO, it is not possible to project only onto the :math:`t_{2g}` subshell because it is not an irreducible representation.
+A redesign of the orthonormalization procedure might happen in the long term, which might allow that.
+
+We strongly suggest using the :py:meth:`.dos_wannier_basis` functionality of the :py:class:`.SumkDFTTools` class (see :download:`calculate_dos_wannier_basis.py <Sr2RuO4/calculate_dos_wannier_basis.py>`) and compare the Wannier-projected orbitals to the original DFT DOS (they should be more or less equal).
+Note that, with SO, there are usually off-diagonal elements of the spectral function, which can also be imaginary.
+The imaginary part can be found in the third column of the files ``DOS_wann_...``.

doc/guide/Sr2RuO4/Sr2RuO4.indmftpr

@@ -0,0 +1,17 @@
+4 ! Nsort
+2 1 2 2 ! Multiplicities
+3 ! lmax
+complex ! Sr
+0 0 0 0
+0 0 0 0
+cubic ! Ru
+0 0 2 0 ! include d-shell as correlated
+0 0 0 0 ! there are no irreps with SO
+1 ! SO-flag
+complex ! O1
+0 0 0 0
+0 0 0 0
+complex ! O2
+0 0 0 0
+0 0 0 0
+-0.7 1.4 ! energy window (Ry)

doc/guide/Sr2RuO4/Sr2RuO4.struct

@@ -0,0 +1,96 @@
+Sr2RuO4                                  s-o calc. M||  0.00  0.00  1.00
+B                            4 39_I
+             RELA
+  7.300012  7.300012 24.044875 90.000000 90.000000 90.000000
+ATOM  -1: X=0.00000000 Y=0.00000000 Z=0.35240000
+          MULT= 2          ISPLIT=-2
+      -1: X=0.00000000 Y=0.00000000 Z=0.64760000
+Sr2+       NPT=  781  R0=.000010000 RMT=   2.26000   Z:  38.00000
+LOCAL ROT MATRIX:    1.0000000 0.0000000 0.0000000
+                     0.0000000 1.0000000 0.0000000
+                     0.0000000 0.0000000 1.0000000
+ATOM  -2: X=0.00000000 Y=0.00000000 Z=0.00000000
+          MULT= 1          ISPLIT=-2
+Ru4+       NPT=  781  R0=.000010000 RMT=   1.95000   Z:  44.00000
+LOCAL ROT MATRIX:    1.0000000 0.0000000 0.0000000
+                     0.0000000 1.0000000 0.0000000
+                     0.0000000 0.0000000 1.0000000
+ATOM  -3: X=0.50000000 Y=0.00000000 Z=0.00000000
+          MULT= 2          ISPLIT= 8
+      -3: X=0.00000000 Y=0.50000000 Z=0.00000000
+O 2-       NPT=  781  R0=.000100000 RMT=   1.68000   Z:   8.00000
+LOCAL ROT MATRIX:    1.0000000 0.0000000 0.0000000
+                     0.0000000 1.0000000 0.0000000
+                     0.0000000 0.0000000 1.0000000
+ATOM  -4: X=0.00000000 Y=0.00000000 Z=0.16350000
+          MULT= 2          ISPLIT=-2
+      -4: X=0.00000000 Y=0.00000000 Z=0.83650000
+O 2-       NPT=  781  R0=.000100000 RMT=   1.68000   Z:   8.00000
+LOCAL ROT MATRIX:    1.0000000 0.0000000 0.0000000
+                     0.0000000 1.0000000 0.0000000
+                     0.0000000 0.0000000 1.0000000
+  16      NUMBER OF SYMMETRY OPERATIONS
+ 0-1 0 0.00000000
+ 1 0 0 0.00000000
+ 0 0-1 0.00000000
+       1   A   2 so. oper.  type  orig. index
+-1 0 0 0.00000000
+ 0-1 0 0.00000000
+ 0 0-1 0.00000000
+       2   A   3
+ 1 0 0 0.00000000
+ 0 1 0 0.00000000
+ 0 0-1 0.00000000
+       3   A   6
+ 0-1 0 0.00000000
+ 1 0 0 0.00000000
+ 0 0 1 0.00000000
+       4   A   8
+ 0 1 0 0.00000000
+-1 0 0 0.00000000
+ 0 0-1 0.00000000
+       5   A   9
+-1 0 0 0.00000000
+ 0-1 0 0.00000000
+ 0 0 1 0.00000000
+       6   A  11
+ 1 0 0 0.00000000
+ 0 1 0 0.00000000
+ 0 0 1 0.00000000
+       7   A  14
+ 0 1 0 0.00000000
+-1 0 0 0.00000000
+ 0 0 1 0.00000000
+       8   A  15
+ 1 0 0 0.00000000
+ 0-1 0 0.00000000
+ 0 0-1 0.00000000
+       9   B   1
+ 0 1 0 0.00000000
+ 1 0 0 0.00000000
+ 0 0-1 0.00000000
+      10   B   4
+ 0-1 0 0.00000000
+-1 0 0 0.00000000
+ 0 0-1 0.00000000
+      11   B   5
+ 1 0 0 0.00000000
+ 0-1 0 0.00000000
+ 0 0 1 0.00000000
+      12   B   7
+-1 0 0 0.00000000
+ 0 1 0 0.00000000
+ 0 0-1 0.00000000
+      13   B  10
+ 0 1 0 0.00000000
+ 1 0 0 0.00000000
+ 0 0 1 0.00000000
+      14   B  12
+ 0-1 0 0.00000000
+-1 0 0 0.00000000
+ 0 0 1 0.00000000
+      15   B  13
+-1 0 0 0.00000000
+ 0 1 0 0.00000000
+ 0 0 1 0.00000000
+      16   B  16

doc/guide/Sr2RuO4/calculate_dos_wannier_basis.py

@@ -0,0 +1,15 @@
+from triqs_dft_tools.converters.wien2k_converter import Wien2kConverter
+from triqs_dft_tools import SumkDFTTools
+
+filename = 'Sr2RuO4'
+
+conv = Wien2kConverter(filename = filename,hdf_filename=filename+'.h5')
+conv.convert_dft_input()
+
+SK = SumkDFTTools(filename+'.h5')
+mesh = (-10.0,10.0,500)
+SK.dos_wannier_basis(broadening=(mesh[1]-mesh[0])/float(mesh[2]),
+                     mesh=mesh,
+                     save_to_file=True,
+                     with_Sigma=False,
+                     with_dc=False)

doc/reference/block_structure.rst

@@ -19,4 +19,3 @@ The block structure can also be written to and read from HDF files.
 .. autoclass:: triqs_dft_tools.block_structure.BlockStructure
    :members:
    :show-inheritance:
-   

python/block_structure.py

@@ -1,9 +1,36 @@
+
+##########################################################################
+#
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+#
+# Copyright (C) 2018 by G. J. Kraberger
+# Copyright (C) 2018 by Simons Foundation
+# Authors: G. J. Kraberger, O. Parcollet
+#
+# TRIQS is free software: you can redistribute it and/or modify it under the
+# terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
+#
+# TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+# details.
+#
+# You should have received a copy of the GNU General Public License along with
+# TRIQS. If not, see <http://www.gnu.org/licenses/>.
+#
+##########################################################################
+
+
 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
+from collections import defaultdict
+
 
 class BlockStructure(object):
     """ Contains information about the Green function structure.
@@ -36,19 +63,268 @@ class BlockStructure(object):
 
         maps from the solver block to the sumk block
         for *inequivalent* correlated shell ish
+    deg_shells : list of lists of lists OR list of lists of dicts
+        In the simple format, ``deg_shells[ish][grp]`` is a list of
+        block names; ``ish`` is the index of the inequivalent correlated shell,
+        ``grp`` is the index of the group of symmetry-related blocks.
+        When the name of two blocks is in the same group, that means that
+        these two blocks are the same by symmetry.
+
+        In the more general format, ``deg_shells[ish][grp]`` is a
+        dictionary whose keys are the block names that are part of the
+        group. The values of the dict for each key are tuples ``(v, conj)``,
+        where ``v`` is a transformation matrix with the same matrix dimensions
+        as the block and ``conj`` is a bool (whether or not to conjugate).
+        Two blocks with ``v_1, conj_1`` and ``v_2, conj_2`` being in the same
+        symmetry group means that
+
+        .. math::
+
+            C_1(v_1^\dagger G_1 v_1) = C_2(v_2^\dagger G_2 v_2),
+
+        where the :math:`G_i` are the Green's functions of the block,
+        and the functions :math:`C_i` conjugate their argument if the bool
+        ``conj_i`` is ``True``.
+    corr_to_inequiv : list
+        a list where, for each correlated shell, the index of the corresponding
+        inequivalent correlated shell is given
+    transformation : list of numpy.array or list of dict
+        a list with entries for each ``ish`` giving transformation matrices
+        that are used on the Green's function in ``sumk`` space when before
+        converting to the ``solver`` space
+        Up to the change in block structure,
+
+        .. math::
+
+            G_{solver} = T G_{sumk} T^\dagger
+
+        if :math:`T` is the ``transformation`` of that particular shell.
+
+        Note that for each shell this can either be a numpy array which
+        applies to all blocks or a dict with a transformation matrix for
+        each block.
     """
-    def __init__(self,gf_struct_sumk=None,
+
+    def __init__(self, gf_struct_sumk=None,
                  gf_struct_solver=None,
                  solver_to_sumk=None,
                  sumk_to_solver=None,
                  solver_to_sumk_block=None,
-                      deg_shells=None):
+                 deg_shells=None,
+                 corr_to_inequiv = None,
+                 transformation=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
+        self.corr_to_inequiv = corr_to_inequiv
+        self.transformation = transformation
+
+    @property
+    def gf_struct_solver_list(self):
+        """ The structure of the solver Green's function
+
+        This is returned as a
+        list (for each shell)
+        of lists (for each block)
+        of tuples (block_name, block_indices).
+
+        That is,
+        ``gf_struct_solver_list[ish][b][0]``
+        is the name of the block number ``b`` of shell ``ish``, and
+        ``gf_struct_solver_list[ish][b][1]``
+        is a list of its indices.
+
+        The list for each shell is sorted alphabetically by block name.
+        """
+        if self.gf_struct_solver is None:
+            return None
+        # we sort by block name in order to get a reproducible result
+        return [sorted([(k, v) for k, v in gfs.iteritems()], key=lambda x: x[0])
+                for gfs in self.gf_struct_solver]
+
+    @property
+    def gf_struct_sumk_list(self):
+        """ The structure of the sumk Green's function
+
+        This is returned as a
+        list (for each shell)
+        of lists (for each block)
+        of tuples (block_name, block_indices)
+
+        That is,
+        ``gf_struct_sumk_list[ish][b][0]``
+        is the name of the block number ``b`` of shell ``ish``, and
+        ``gf_struct_sumk_list[ish][b][1]``
+        is a list of its indices.
+        """
+        return self.gf_struct_sumk
+
+    @property
+    def gf_struct_solver_dict(self):
+        """ The structure of the solver Green's function
+
+        This is returned as a
+        list (for each shell)
+        of dictionaries.
+
+        That is,
+        ``gf_struct_solver_dict[ish][bname]``
+        is a list of the indices of block ``bname`` of shell ``ish``.
+        """
+        return self.gf_struct_solver
+
+    @property
+    def gf_struct_sumk_dict(self):
+        """ The structure of the sumk Green's function
+
+        This is returned as a
+        list (for each shell)
+        of dictionaries.
+
+        That is,
+        ``gf_struct_sumk_dict[ish][bname]``
+        is a list of the indices of block ``bname`` of shell ``ish``.
+        """
+        if self.gf_struct_sumk is None:
+            return None
+        return [{block: indices for block, indices in gfs}
+                for gfs in self.gf_struct_sumk]
+
+    @property
+    def inequiv_to_corr(self):
+        """ A list mapping an inequivalent correlated shell to a correlated shell
+        """
+
+        if self.corr_to_inequiv is None:
+            return None
+        N_solver = len(np.unique(self.corr_to_inequiv))
+        if self.gf_struct_solver is not None:
+            assert N_solver == len(self.gf_struct_solver)
+        assert sorted(np.unique(self.corr_to_inequiv)) == range(N_solver),\
+            "an inequivalent shell is missing in corr_to_inequiv"
+        return [self.corr_to_inequiv.index(icrsh)
+                for icrsh in range(N_solver)]
+
+    @inequiv_to_corr.setter
+    def inequiv_to_corr(self, value):
+        # a check for backward compatibility
+        if value is None:
+            return
+        assert self.inequiv_to_corr == value, "Trying to set incompatible inequiv_to_corr"
+
+    @property
+    def sumk_to_solver_block(self):
+        if self.inequiv_to_corr is None:
+            return None
+        ret = []
+        for ish, icrsh in enumerate(self.inequiv_to_corr):
+            d = defaultdict(list)
+            for block_solver, block_sumk in self.solver_to_sumk_block[ish].iteritems():
+                d[block_sumk].append(block_solver)
+            ret.append(d)
+        return ret
+
+    @property
+    def effective_transformation_sumk(self):
+        """ Return the effective transformation matrix
+
+        A list of dicts, one for every correlated shell. In the dict,
+        there is a transformation matrix (as numpy array) for each
+        block in sumk space, that is used to transform the block.
+        """
+        trans = copy.deepcopy(self.transformation)
+        if self.gf_struct_sumk is None:
+            raise Exception('gf_struct_sumk not set.')
+        if self.gf_struct_solver is None:
+            raise Exception('gf_struct_solver not set.')
+
+        if trans is None:
+            trans = [{block: np.eye(len(indices)) for block, indices in gfs}
+                     for gfs in self.gf_struct_sumk]
+
+        assert isinstance(trans, list),\
+            "transformation has to be a list"
+
+        assert len(trans) == len(self.gf_struct_sumk),\
+            "give one transformation per correlated shell"
+
+        for icrsh in range(len(trans)):
+            ish = self.corr_to_inequiv[icrsh]
+            if trans[icrsh] is None:
+                trans[icrsh] = {block: np.eye(len(indices))
+                                for block, indices in self.gf_struct_sumk[icrsh]}
+
+            if not isinstance(trans[icrsh], dict):
+                trans[icrsh] = {block: copy.deepcopy(trans[icrsh])
+                                for block, indices in self.gf_struct_sumk[icrsh]}
+
+            assert trans[icrsh].keys() == self.gf_struct_sumk_dict[icrsh].keys(),\
+                "wrong block names used in transformation (icrsh = {})".format(icrsh)
+
+            for block in trans[icrsh]:
+                assert trans[icrsh][block].shape[0] == trans[icrsh][block].shape[1],\
+                    "Transformation has to be quadratic; throwing away orbitals can be achieved on the level of the mapping. (icrsh = {}, block = {})".format(icrsh, block)
+
+                assert trans[icrsh][block].shape[0] == len(self.gf_struct_sumk_dict[icrsh][block]),\
+                    "Transformation block shape does not match gf_struct_sumk. (icrsh = {}, block = {})".format(icrsh, block)
+
+                # zero out all the lines of the transformation that are
+                # not included in gf_struct_solver
+                for iorb, norb in enumerate(self.gf_struct_sumk_dict[icrsh][block]):
+                    if self.sumk_to_solver[ish][(block, norb)][0] is None:
+                        trans[icrsh][block][iorb, :] = 0.0
+        return trans
+
+    @property
+    def effective_transformation_solver(self):
+        """ Return the effective transformation matrix
+
+        A list of dicts, one for every inequivalent correlated shell.
+        In the dict, there is a transformation matrix (as numpy array)
+        for each block in solver space, that is used to transform from
+        the sumk block (see :py:meth:`.solver_to_sumk_block`) to the
+        solver block.
+
+
+        For a solver block ``b`` for inequivalent correlated shell ``ish``,
+        the corresponding block of the solver Green's function is::
+
+            # the effective transformation matrix for the block
+            T = block_structure.effective_transformation_solver[ish][b]
+            # the index of the correlated shell
+            icrsh = block_structure.inequiv_to_corr[ish]
+            # the name of the corresponding sumk block
+            block_sumk = block_structure.solver_to_sumk_block[icrsh][b]
+            # transform the Green's function
+            G_solver[ish][b].from_L_G_R(T, G_sumk[icrsh][block_sumk], T.conjugate().transpose())
+
+        The functionality of that code block is implemented in
+        :py:meth:`.convert_gf` (i.e., you don't need to use this directly).
+        """
+
+        eff_trans_sumk = self.effective_transformation_sumk
+
+        ets = []
+        for ish in range(len(self.gf_struct_solver)):
+            icrsh = self.inequiv_to_corr[ish]
+            ets.append(dict())
+            for block in self.gf_struct_solver[ish]:
+                block_sumk = self.solver_to_sumk_block[ish][block]
+                T = eff_trans_sumk[icrsh][block_sumk]
+                ets[ish][block] = np.zeros((len(self.gf_struct_solver[ish][block]),
+                                            len(T)),
+                                           dtype=T.dtype)
+                for i in self.gf_struct_solver[ish][block]:
+                    i_sumk = self.solver_to_sumk[ish][block, i]
+                    assert i_sumk[0] == block_sumk,\
+                        "Wrong block in solver_to_sumk"
+                    i_sumk = i_sumk[1]
+                    ets[ish][block][i, :] = T[i_sumk, :]
+        return ets
+
 
     @classmethod
     def full_structure(cls,gf_struct,corr_to_inequiv):
@@ -103,9 +379,10 @@ class BlockStructure(object):
                 solver_to_sumk = copy.deepcopy(solver_to_sumk),
                 sumk_to_solver = solver_to_sumk,
                 solver_to_sumk_block = s2sblock,
-                deg_shells = [[] for ish in range(len(gf_struct))])
+                deg_shells = [[] for ish in range(len(gf_struct))],
+                corr_to_inequiv = corr_to_inequiv)
 
-    def pick_gf_struct_solver(self,new_gf_struct):
+    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
@@ -143,28 +420,57 @@ class BlockStructure(object):
             gf_struct = new_gf_struct[ish]
 
             # create new solver_to_sumk
-            so2su={}
+            so2su = {}
             so2su_block = {}
-            for blk,idxs in gf_struct.items():
+            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[(blk, i)] = self.solver_to_sumk[ish][(blk, idxs[i])]
-                    so2su_block[blk]=so2su[(blk,i)][0]
+                    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():
+            for k, v in self.sumk_to_solver[ish].items():
-                blk,ind=v
+                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)
+                    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)
+            # adapt deg_shells
+
+            self.adapt_deg_shells(gf_struct, ish)
+
+
             # reindexing gf_struct so that it starts with 0
             for k in gf_struct:
-                gf_struct[k]=range(len(gf_struct[k]))
+                gf_struct[k] = range(len(gf_struct[k]))
-            self.gf_struct_solver[ish]=gf_struct
+            self.gf_struct_solver[ish] = gf_struct
 
-    def pick_gf_struct_sumk(self,new_gf_struct):
+    def adapt_deg_shells(self, gf_struct, ish=0):
+        """ Adapts the deg_shells to a new gf_struct
+            Internally called when using pick_gf_struct and map_gf_struct
+        """
+        if self.deg_shells is not None:
+            for degsh in self.deg_shells[ish]:
+                if isinstance(degsh, dict):
+                    for key in degsh.keys():
+                        if not key in gf_struct:
+                            del degsh[key]
+                            continue
+                        if gf_struct[key] != self.gf_struct_solver[ish][key]:
+                            v, C = degsh[key]
+                            degsh[key][0] = \
+                                v[gf_struct[key], :][:, gf_struct[key]]
+                            warn(
+                                'Removing shells from degenerate shell {}. Cannot guarantee that they continue to be equivalent.')
+                else: # degshell is list
+                    degsh1 = copy.deepcopy(degsh) # in order to not remove a key while iterating
+                    for key in degsh1:
+                        if not key in gf_struct:
+                            warn('Removing shells from degenerate shell {}.')
+                            degsh.remove(key)
+
+    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
@@ -199,25 +505,50 @@ class BlockStructure(object):
             However, the indices are not according to the solver Gf
             but the sumk Gf.
         """
+        eff_trans_sumk = self.effective_transformation_sumk
 
+        assert len(eff_trans_sumk) == len(new_gf_struct),\
+            "new_gf_struct has the wrong length"
+
+        new_gf_struct_transformed = copy.deepcopy(new_gf_struct)
+
+        # when there is a transformation matrix, this first zeroes out
+        # the corresponding rows of (a copy of) T and then applies
+        # pick_gf_struct_solver for all lines of  T that have at least
+        # one non-zero entry
+
+        for icrsh in range(len(new_gf_struct)):
+            for block, indices in self.gf_struct_sumk[icrsh]:
+                if not block in new_gf_struct[icrsh]:
+                    #del new_gf_struct_transformed[block] # this MUST be wrong, as new_gf_struct_transformed needs to have a integer index for icrsh... # error when index is not kept at all
+                    continue
+                T = eff_trans_sumk[icrsh][block]
+                for index in indices:
+                    if not index in new_gf_struct[icrsh][block]:
+                        T[:, index] = 0.0
+                new_indices = []
+                for index in indices:
+                    if np.any(np.abs(T[index, :]) > 1.e-15):
+                        new_indices.append(index)
+                new_gf_struct_transformed[icrsh][block] = new_indices
 
         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)):
+        for icrsh in range(len(new_gf_struct_transformed)):
+            ish = self.corr_to_inequiv[icrsh]
             gfs.append({})
-            for block in new_gf_struct[ish].keys():
+            for block in new_gf_struct_transformed[icrsh].keys():
-                for ind in new_gf_struct[ish][block]:
+                for ind in new_gf_struct_transformed[icrsh][block]:
-                    ind_sol = self.sumk_to_solver[ish][(block,ind)]
+                    ind_sol = self.sumk_to_solver[ish][(block, ind)]
-                    if not ind_sol[0] in gfs[ish]:
+                    if not ind_sol[0] in gfs[icrsh]:
-                        gfs[ish][ind_sol[0]]=[]
+                        gfs[icrsh][ind_sol[0]] = []
-                    gfs[ish][ind_sol[0]].append(ind_sol[1])
+                    gfs[icrsh][ind_sol[0]].append(ind_sol[1])
         self.pick_gf_struct_solver(gfs)
 
-
+    def map_gf_struct_solver(self, mapping):
-    def map_gf_struct_solver(self,mapping):
+        r""" Map the Green function structure from one struct to another.
-        """ Map the Green function structure from one struct to another.
 
         Parameters
         ----------
@@ -225,38 +556,61 @@ class BlockStructure(object):
             the dict consists of elements
             (from_block,from_index) : (to_block,to_index)
             that maps from one structure to the other
+            (one for each shell; use a mapping ``None`` for a shell
+            you want to leave unchanged)
+
+        Examples
+        --------
+
+        Consider a `gf_struct_solver` consisting of two :math:`1 \times 1`
+        blocks, `block_1` and `block_2`. Say you want to have a new block
+        structure where you merge them into one block because you want to
+        introduce an off-diagonal element. You could perform the mapping
+        via::
+
+            map_gf_struct_solver([{('block_1',0) : ('block', 0)
+                                   ('block_2',0) : ('block', 1)}])
         """
 
         for ish in range(len(mapping)):
+            if mapping[ish] is None:
+                continue
             gf_struct = {}
             so2su = {}
             su2so = {}
             so2su_block = {}
-            for frm,to in mapping[ish].iteritems():
+            for frm, to in mapping[ish].iteritems():
                 if not to[0] in gf_struct:
-                    gf_struct[to[0]]=[]
+                    gf_struct[to[0]] = []
                 gf_struct[to[0]].append(to[1])
 
-                so2su[to]=self.solver_to_sumk[ish][frm]
+                so2su[to] = self.solver_to_sumk[ish][frm]
-                su2so[self.solver_to_sumk[ish][frm]]=to
+                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]]))
+                                to[0], so2su_block[to[0]], self.solver_to_sumk_block[ish][frm[0]]))
                 else:
-                    so2su_block[to[0]]=\
+                    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)
+                    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):
+            self.adapt_deg_shells(gf_struct, ish)
-        """ Create a zero BlockGf having the gf_struct_solver structure.
+
+            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
+            self.deg_shells[ish] = []
+
+    def create_gf(self, ish=0, gf_function=GfImFreq, space='solver', **kwargs):
+        """ Create a zero BlockGf having the correct structure.
+
+        For ``space='solver'``, the structure is according to
+        ``gf_struct_solver``, else according to ``gf_struct_sumk``.
 
         When using GfImFreq as gf_function, typically you have to
         supply beta as keyword argument.
@@ -265,24 +619,158 @@ class BlockStructure(object):
         ----------
         ish : int
             shell index
+            If ``space='solver'``, the index of the of the inequivalent correlated shell,
+            if ``space='sumk'``, the index of the correlated shell
         gf_function : constructor
             function used to construct the Gf objects constituting the
             individual blocks; default: GfImFreq
+        space : 'solver' or 'sumk'
+            which space the structure should correspond to
         **kwargs :
             options passed on to the Gf constructor for the individual
             blocks
         """
 
-        names = self.gf_struct_solver[ish].keys()
+        return self._create_gf_or_matrix(ish, gf_function, BlockGf, space, **kwargs)
-        blocks=[]
+
+    def create_matrix(self, ish=0, space='solver', dtype=np.complex_):
+        """ Create a zero matrix having the correct structure.
+
+        For ``space='solver'``, the structure is according to
+        ``gf_struct_solver``, else according to ``gf_struct_sumk``.
+
+        Parameters
+        ----------
+        ish : int
+            shell index
+            If ``space='solver'``, the index of the of the inequivalent correlated shell,
+            if ``space='sumk'``, the index of the correlated shell
+        space : 'solver' or 'sumk'
+            which space the structure should correspond to
+        """
+
+        def gf_function(indices):
+            return np.zeros((len(indices), len(indices)), dtype=dtype)
+
+        def block_function(name_list, block_list):
+            d = dict()
+            for i in range(len(name_list)):
+                d[name_list[i]] = block_list[i]
+            return d
+
+        return self._create_gf_or_matrix(ish, gf_function, block_function, space)
+
+    def _create_gf_or_matrix(self, ish=0, gf_function=GfImFreq, block_function=BlockGf, space='solver', **kwargs):
+        if space == 'solver':
+            gf_struct = self.gf_struct_solver
+        elif space == 'sumk':
+            gf_struct = self.gf_struct_sumk_dict
+        else:
+            raise Exception(
+                "Argument space has to be either 'solver' or 'sumk'.")
+
+        names = gf_struct[ish].keys()
+        blocks = []
         for n in names:
-            G = gf_function(indices=self.gf_struct_solver[ish][n],**kwargs)
+            G = gf_function(indices=gf_struct[ish][n], **kwargs)
             blocks.append(G)
-        G = BlockGf(name_list = names, block_list = blocks)
+        G = block_function(name_list=names, block_list=blocks)
         return G
 
+    def check_gf(self, G, ish=None, space='solver'):
+        """ check whether the Green's function G has the right structure
 
-    def convert_gf(self,G,G_struct,ish=0,show_warnings=True,**kwargs):
+        This throws an error if the structure of G is not the same
+        as ``gf_struct_solver`` (for ``space=solver``) or
+        ``gf_struct_sumk`` (for ``space=sumk``)..
+
+        Parameters
+        ----------
+        G : BlockGf or list of BlockGf
+            Green's function to check
+            if it is a list, there should be as many entries as there
+            are shells, and the check is performed for all shells (unless
+            ish is given).
+        ish : int
+            shell index
+            default: 0 if G is just one Green's function is given,
+            check all if list of Green's functions is given
+        space : 'solver' or 'sumk'
+            which space the structure should correspond to
+        """
+
+        return self._check_gf_or_matrix(G, ish, space)
+
+    def check_matrix(self, G, ish=None, space='solver'):
+        """ check whether the matrix G has the right structure
+
+        This throws an error if the structure of G is not the same
+        as ``gf_struct_solver`` (for ``space=solver``) or
+        ``gf_struct_sumk`` (for ``space=sumk``)..
+
+        Parameters
+        ----------
+        G : dict of matrices or list of dict of matrices
+            matrix to check
+            if it is a list, there should be as many entries as there
+            are shells, and the check is performed for all shells (unless
+            ish is given).
+        ish : int
+            shell index
+            default: 0 if G is just one matrix is given,
+            check all if list of dicts is given
+        space : 'solver' or 'sumk'
+            which space the structure should correspond to
+        """
+
+        return self._check_gf_or_matrix(G, ish, space)
+
+    def _check_gf_or_matrix(self, G, ish=None, space='solver'):
+        if space == 'solver':
+            gf_struct = self.gf_struct_solver
+        elif space == 'sumk':
+            gf_struct = self.gf_struct_sumk_dict
+        else:
+            raise Exception(
+                "Argument space has to be either 'solver' or 'sumk'.")
+
+        if isinstance(G, list):
+            assert len(G) == len(gf_struct),\
+                "list of G does not have the correct length"
+            if ish is None:
+                ishs = range(len(gf_struct))
+            else:
+                ishs = [ish]
+            for ish in ishs:
+                self.check_gf(G[ish], ish=ish, space=space)
+            return
+
+        if ish is None:
+            ish = 0
+
+        if isinstance(G, BlockGf):
+            for block in gf_struct[ish]:
+                assert block in G.indices,\
+                    "block " + block + " not in G (shell {})".format(ish)
+            for block, gf in G:
+                assert block in gf_struct[ish],\
+                    "block " + block + " not in struct (shell {})".format(ish)
+                assert list(gf.indices) == 2 * [map(str, gf_struct[ish][block])],\
+                    "block " + block + \
+                    " has wrong indices (shell {})".format(ish)
+        else:
+            for block in gf_struct[ish]:
+                assert block in G,\
+                    "block " + block + " not in G (shell {})".format(ish)
+            for block, gf in G.iteritems():
+                assert block in gf_struct[ish],\
+                    "block " + block + " not in struct (shell {})".format(ish)
+                assert range(len(gf)) == 2 * [map(str, gf_struct[ish][block])],\
+                    "block " + block + \
+                    " has wrong indices (shell {})".format(ish)
+
+    def convert_gf(self, G, G_struct=None, ish_from=0, ish_to=None, show_warnings=True,
+                   G_out=None, space_from='solver', space_to='solver', ish=None, **kwargs):
         """ Convert BlockGf from its structure to this structure.
 
         .. warning::
@@ -295,49 +783,194 @@ class BlockStructure(object):
         ----------
         G : BlockGf
             the Gf that should be converted
-        G_struct : GfStructure
+        G_struct : BlockStructure or str
-            the structure of that G
+            the structure of that G or None (then, this structure
-        ish : int
+            is used)
-            shell index
+        ish_from : int
+            shell index of the input structure
+        ish_to : int
+            shell index of the output structure; if None (the default),
+            it is the same as ish_from
         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)
+        G_out : BlockGf
+            the output Green's function (if not given, a new one is
+            created)
+        space_from : 'solver' or 'sumk'
+            whether the Green's function ``G`` corresponds to the
+            solver or sumk structure of ``G_struct``
+        space_to : 'solver' or 'sumk'
+            whether the output Green's function should be according to
+            the solver of sumk structure of this structure
         **kwargs :
             options passed to the constructor for the new Gf
         """
 
+        if ish is not None:
+            warn(
+                'The parameter ish in convert_gf is deprecated. Use ish_from and ish_to instead.')
+            ish_from = ish
+            ish_to = ish
+        return self._convert_gf_or_matrix(G, G_struct, ish_from, ish_to,
+                                          show_warnings, G_out, space_from, space_to, **kwargs)
+
+    def convert_matrix(self, G, G_struct=None, ish_from=0, ish_to=None, show_warnings=True,
+                       G_out=None, space_from='solver', space_to='solver'):
+        """ Convert matrix 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 : dict of numpy array
+            the matrix that should be converted
+        G_struct : BlockStructure or str
+            the structure of that G or None (then, this structure
+            is used)
+        ish_from : int
+            shell index of the input structure
+        ish_to : int
+            shell index of the output structure; if None (the default),
+            it is the same as ish_from
+        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)
+        G_out : dict of numpy array
+            the output numpy array (if not given, a new one is
+            created)
+        space_from : 'solver' or 'sumk'
+            whether the matrix ``G`` corresponds to the
+            solver or sumk structure of ``G_struct``
+        space_to : 'solver' or 'sumk'
+            whether the output matrix should be according to
+            the solver of sumk structure of this structure
+        **kwargs :
+            options passed to the constructor for the new Gf
+        """
+
+        return self._convert_gf_or_matrix(G, G_struct, ish_from, ish_to,
+                                          show_warnings, G_out, space_from, space_to)
+
+    def _convert_gf_or_matrix(self, G, G_struct=None, ish_from=0, ish_to=None, show_warnings=True,
+                              G_out=None, space_from='solver', space_to='solver', **kwargs):
+        if ish_to is None:
+            ish_to = ish_from
+
         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)
+        if G_struct is None:
-        for block in G_struct.gf_struct_solver[ish].keys():
+            G_struct = self
-            for i1 in G_struct.gf_struct_solver[ish][block]:
+
-                for i2 in G_struct.gf_struct_solver[ish][block]:
+        if space_from == 'solver':
-                    i1_sumk = G_struct.solver_to_sumk[ish][(block,i1)]
+            gf_struct_from = G_struct.gf_struct_solver[ish_from]
-                    i2_sumk = G_struct.solver_to_sumk[ish][(block,i2)]
+            eff_trans_from = G_struct.effective_transformation_solver[ish_from]
-                    i1_sol = self.sumk_to_solver[ish][i1_sumk]
+            block_mapping_from = G_struct.sumk_to_solver_block[ish_from]
-                    i2_sol = self.sumk_to_solver[ish][i2_sumk]
+        elif space_from == 'sumk':
-                    if i1_sol[0] is None or i2_sol[0] is None:
+            gf_struct_from = G_struct.gf_struct_sumk_dict[ish_from]
+            eff_trans_from = {block: np.eye(len(indices))
+                              for block, indices in G_struct.gf_struct_sumk[ish_from]}
+            block_mapping_from = {b: [b] for b in gf_struct_from}
+        else:
+            raise Exception(
+                "Argument space_from has to be either 'solver' or 'sumk'.")
+
+        if space_to == 'solver':
+            gf_struct_to = self.gf_struct_solver[ish_to]
+            eff_trans_to = self.effective_transformation_solver[ish_to]
+            block_mapping_to = self.solver_to_sumk_block[ish_to]
+        elif space_to == 'sumk':
+            gf_struct_to = self.gf_struct_sumk_dict[ish_to]
+            eff_trans_to = {block: np.eye(len(indices))
+                            for block, indices in self.gf_struct_sumk_list[ish_to]}
+            block_mapping_to = {b: b for b in gf_struct_to}
+        else:
+            raise Exception(
+                "Argument space_to has to be either 'solver' or 'sumk'.")
+
+        if isinstance(G, BlockGf):
+            # create a Green's function to hold the result
+            if G_out is None:
+                if not 'mesh' in kwargs and not 'beta' in kwargs:
+                    kwargs['mesh'] = G.mesh
+                G_out = self.create_gf(ish=ish_to, space=space_to, **kwargs)
+            else:
+                self.check_gf(G_out, ish=ish_to, space=space_to)
+        elif isinstance(G, dict):
+            if G_out is None:
+                G_out = self.create_matrix(ish=ish_to, space=space_to)
+            else:
+                self.check_matrix(G_out, ish=ish_to, space=space_to)
+        else:
+            raise Exception('G is neither BlockGf nor dict.')
+
+        for block_to in gf_struct_to.keys():
+            if isinstance(G, BlockGf):
+                G_out[block_to].zero()
+            else:
+                G_out[block_to][:] = 0.0
+            block_intermediate = block_mapping_to[block_to]
+            block_from = block_mapping_from[block_intermediate]
+            T_to = eff_trans_to[block_to]
+            g_help = G_out[block_to].copy()
+            for block in block_from:
+                T_from = eff_trans_from[block]
+                if isinstance(G, BlockGf):
+                    g_help.from_L_G_R(np.dot(T_to, np.conjugate(np.transpose(T_from))),
+                                      G[block],
+                                      np.dot(T_from, np.conjugate(np.transpose(T_to))))
+                    G_out[block_to] << G_out[block_to] + g_help
+                else:
+                    g_help = np.dot(np.dot(T_to, np.conjugate(np.transpose(T_from))),
+                                    np.dot(G[block],
+                                           np.dot(T_from, np.conjugate(np.transpose(T_to)))))
+                    G_out[block_to] += g_help
+
         if show_warnings:
-                            if mpi.is_master_node():
+            # we back-transform it
-                                warn(('Element {},{} of block {} of G is not present '+
+            G_back = G_struct._convert_gf_or_matrix(G_out, self, ish_from=ish_to,
-                                    'in the new structure').format(i1,i2,block))
+                                                    ish_to=ish_from,
-                        continue
+                                                    show_warnings=False,  # else we get an endless loop
-                    if i1_sol[0]!=i2_sol[0]:
+                                                    space_from=space_to, space_to=space_from, **kwargs)
-                        if show_warnings and np.max(np.abs(G[block][i1,i2].data)) > warning_threshold:
+            for name, gf in (G_back if isinstance(G, BlockGf) else G_back.iteritems()):
-                            if mpi.is_master_node():
+                if isinstance(G, BlockGf):
-                                warn(('Element {},{} of block {} of G is approximated '+
+                    maxdiff = np.max(np.abs(G_back[name].data - G[name].data),
-                                    'to zero to match the new structure. Max abs value: {}').format(
+                                     axis=0)
-                                        i1,i2,block,np.max(np.abs(G[block][i1,i2].data))))
+                else:
-                        continue
+                    maxdiff = G_back[name] - G[name]
-                    G_new[i1_sol[0]][i1_sol[1],i2_sol[1]] = \
+
-                            G[block][i1,i2]
+                if space_to == 'solver' and self == G_struct: # do comparison in solver (ignore diff. in ignored orbitals)
-        return G_new
+                    tmp = self.create_matrix(space='sumk')
+                    tmp[name] = maxdiff
+                    maxdiff = G_struct._convert_gf_or_matrix(tmp, self, ish_from=ish_from,
+                                                    ish_to=ish_to,
+                                                    show_warnings=False,
+                                                    space_from=space_from, space_to=space_to, **kwargs)
+
+                    for block in maxdiff:
+                        maxdiff_b = maxdiff[block]
+                        if np.any(maxdiff_b > warning_threshold):
+                            warn('Block {} maximum difference:\n'.format(name) + str(maxdiff))
+
+
+                elif np.any(maxdiff > warning_threshold):
+                    warn('Block {} maximum difference:\n'.format(name)
+                         + str(maxdiff))
+
+        return G_out
 
     def approximate_as_diagonal(self):
         """ Create a structure for a GF with zero off-diagonal elements.
@@ -393,7 +1026,7 @@ class BlockStructure(object):
 
         for prop in [ "gf_struct_sumk", "gf_struct_solver",
                 "solver_to_sumk", "sumk_to_solver", "solver_to_sumk_block",
-                "deg_shells"]:
+                "deg_shells","transformation", "corr_to_inequiv"]:
             if not compare(getattr(self,prop),getattr(other,prop)):
                 return False
         return True
@@ -406,8 +1039,14 @@ class BlockStructure(object):
 
         ret = {}
         for element in [ "gf_struct_sumk", "gf_struct_solver",
-                         "solver_to_sumk_block","deg_shells"]:
+                         "solver_to_sumk_block","deg_shells",
+                         "transformation", "corr_to_inequiv"]:
             ret[element] = getattr(self,element)
+            if ret[element] is None:
+                ret[element] = 'None'
+
+        if ret["transformation"] is None:
+            ret["transformation"] = "None"
 
         def construct_mapping(mapping):
             d = []
@@ -434,13 +1073,18 @@ class BlockStructure(object):
                     d[ish][literal_eval(k)] = literal_eval(v)
             return d
 
+        for elem in D:
+            if D[elem]=="None":
+                D[elem] = None
+
         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+= "corr_to_inequiv "+str(self.corr_to_inequiv)+'\n'
+        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']:
@@ -465,6 +1109,8 @@ class BlockStructure(object):
                 else:
                     for key in self.deg_shells[ish][l]:
                         s+='   '+key+'\n'
+        s += "transformation\n"
+        s += str(self.transformation)
         return s
 
 from pytriqs.archive.hdf_archive_schemes import register_class

python/sumk_dft.py

@@ -3,6 +3,7 @@
 #
 # TRIQS: a Toolbox for Research in Interacting Quantum Systems
 #
+# Copyright (C) 2018 by G. J. Kraberger
 # Copyright (C) 2011 by M. Aichhorn, L. Pourovskii, V. Vildosola
 #
 # TRIQS is free software: you can redistribute it and/or modify it under the
@@ -94,6 +95,8 @@ class SumkDFT(object):
             self.misc_data = misc_data
             self.h_field = h_field
 
+            self.block_structure = BlockStructure()
+
             # Read input from HDF:
             things_to_read = ['energy_unit', 'n_k', 'k_dep_projection', 'SP', 'SO', 'charge_below', 'density_required',
                               'symm_op', 'n_shells', 'shells', 'n_corr_shells', 'corr_shells', 'use_rotations', 'rot_mat',
@@ -119,8 +122,6 @@ class SumkDFT(object):
                     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
@@ -183,7 +184,7 @@ class SumkDFT(object):
         # initialise variables on all nodes to ensure mpi broadcast works at
         # the end
         for it in things_to_read:
-            setattr(self, it, 0)
+            setattr(self, it, None)
         subgroup_present = 0
 
         if mpi.is_master_node():
@@ -571,63 +572,98 @@ class SumkDFT(object):
 
         return G_latt
 
-    def set_Sigma(self, Sigma_imp):
+    def set_Sigma(self, Sigma_imp, transform_to_sumk_blocks=True):
-        self.put_Sigma(Sigma_imp)
+        self.put_Sigma(Sigma_imp, transform_to_sumk_blocks)
 
-    def put_Sigma(self, Sigma_imp):
+    def put_Sigma(self, Sigma_imp, transform_to_sumk_blocks=True):
         r"""
-        Inserts the impurity self-energies into the sumk_dft class.
+        Insert the impurity self-energies into the sumk_dft class.
+
+        Parameters
+        ----------
+        Sigma_imp : list of BlockGf (Green's function) objects
+            List containing impurity self-energy for all (inequivalent) correlated shells.
+            Self-energies for equivalent shells are then automatically set by this function.
+            The self-energies can be of the real or imaginary-frequency type.
+        transform_to_sumk_blocks : bool, optional
+            If True (default), the input Sigma_imp will be transformed to the block structure ``gf_struct_sumk``,
+            else it has to be given in ``gf_struct_sumk``.
+        """
+
+        if transform_to_sumk_blocks:
+            Sigma_imp = self.transform_to_sumk_blocks(Sigma_imp)
+
+        assert isinstance(Sigma_imp, list),\
+            "put_Sigma: Sigma_imp has to be a list of Sigmas for the correlated shells, even if it is of length 1!"
+        assert len(Sigma_imp) == self.n_corr_shells,\
+            "put_Sigma: give exactly one Sigma for each corr. shell!"
+
+        if all((isinstance(gf, Gf) and isinstance(gf.mesh, MeshImFreq)) for bname, gf in Sigma_imp[0]):
+            # Imaginary frequency Sigma:
+            self.Sigma_imp_iw = [self.block_structure.create_gf(ish=icrsh, mesh=Sigma_imp[icrsh].mesh, space='sumk')
+                                 for icrsh in range(self.n_corr_shells)]
+            SK_Sigma_imp = self.Sigma_imp_iw
+        elif all(isinstance(gf, Gf) and isinstance(gf.mesh, MeshReFreq) for bname, gf in Sigma_imp[0]):
+            # Real frequency Sigma:
+            self.Sigma_imp_w = [self.block_structure.create_gf(ish=icrsh, mesh=Sigma_imp[icrsh].mesh, gf_function=GfReFreq, space='sumk')
+                                for icrsh in range(self.n_corr_shells)]
+            SK_Sigma_imp = self.Sigma_imp_w
+        else:
+            raise ValueError, "put_Sigma: This type of Sigma is not handled, give either BlockGf of GfReFreq or GfImFreq."
+
+        # rotation from local to global coordinate system:
+        for icrsh in range(self.n_corr_shells):
+            for bname, gf in SK_Sigma_imp[icrsh]:
+                if self.use_rotations:
+                    gf << self.rotloc(icrsh,
+                                      Sigma_imp[icrsh][bname],
+                                      direction='toGlobal')
+                else:
+                    gf << Sigma_imp[icrsh][bname]
+
+    def transform_to_sumk_blocks(self, Sigma_imp, Sigma_out=None):
+        r""" transform Sigma from solver to sumk space
 
         Parameters
         ----------
         Sigma_imp : list of BlockGf (Green's function) objects
             List containing impurity self-energy for all inequivalent correlated shells.
-                    Self-energies for equivalent shells are then automatically set by this function.
             The self-energies can be of the real or imaginary-frequency type.
+        Sigma_out : list of BlockGf
+            list of one BlockGf per correlated shell with the block structure
+            according to ``gf_struct_sumk``; if None, it will be created
         """
 
-        assert isinstance(
+        assert isinstance(Sigma_imp, list),\
-            Sigma_imp, list), "put_Sigma: Sigma_imp has to be a list of Sigmas for the correlated shells, even if it is of length 1!"
+            "transform_to_sumk_blocks: Sigma_imp has to be a list of Sigmas for the inequivalent correlated shells, even if it is of length 1!"
-        assert len(
+        assert len(Sigma_imp) == self.n_inequiv_shells,\
-            Sigma_imp) == self.n_inequiv_shells, "put_Sigma: give exactly one Sigma for each inequivalent corr. shell!"
+            "transform_to_sumk_blocks: give exactly one Sigma for each inequivalent corr. shell!"
 
-        # init self.Sigma_imp_(i)w:
+        if Sigma_out is None:
-        if all( (isinstance(gf, Gf) and isinstance (gf.mesh, MeshImFreq)) for bname, gf in Sigma_imp[0]):
+            Sigma_out = [self.block_structure.create_gf(ish=icrsh, mesh=Sigma_imp[self.corr_to_inequiv[icrsh]].mesh, space='sumk')
-            # Imaginary frequency Sigma:
-            self.Sigma_imp_iw = [BlockGf(name_block_generator=[(block, GfImFreq(indices=inner, mesh=Sigma_imp[0].mesh))
-                                                               for block, inner in self.gf_struct_sumk[icrsh]], make_copies=False)
                          for icrsh in range(self.n_corr_shells)]
-            SK_Sigma_imp = self.Sigma_imp_iw
-        elif all( isinstance(gf, Gf) and isinstance (gf.mesh, MeshReFreq) for bname, gf in Sigma_imp[0]):
-            # Real frequency Sigma:
-            self.Sigma_imp_w = [BlockGf(name_block_generator=[(block, GfReFreq(indices=inner, mesh=Sigma_imp[0].mesh))
-                                                              for block, inner in self.gf_struct_sumk[icrsh]], make_copies=False)
-                                for icrsh in range(self.n_corr_shells)]
-            SK_Sigma_imp = self.Sigma_imp_w
         else:
-            raise ValueError, "put_Sigma: This type of Sigma is not handled."
+            for icrsh in range(self.n_corr_shells):
+                self.block_structure.check_gf(Sigma_out,
+                                              ish=icrsh,
+                                              space='sumk')
 
         # transform the CTQMC blocks to the full matrix:
         for icrsh in range(self.n_corr_shells):
             # ish is the index of the inequivalent shell corresponding to icrsh
             ish = self.corr_to_inequiv[icrsh]
-            for block, inner in self.gf_struct_solver[ish].iteritems():
+            self.block_structure.convert_gf(
-                for ind1 in inner:
+                G=Sigma_imp[ish],
-                    for ind2 in inner:
+                G_struct=None,
-                        block_sumk, ind1_sumk = self.solver_to_sumk[
+                space_from='solver',
-                            ish][(block, ind1)]
+                space_to='sumk',
-                        block_sumk, ind2_sumk = self.solver_to_sumk[
+                ish_from=ish,
-                            ish][(block, ind2)]
+                ish_to=icrsh,
-                        SK_Sigma_imp[icrsh][block_sumk][
+                G_out=Sigma_out[icrsh])
-                            ind1_sumk, ind2_sumk] << Sigma_imp[ish][block][ind1, ind2]
+        return Sigma_out
 
-        # rotation from local to global coordinate system:
+    def extract_G_loc(self, mu=None, iw_or_w='iw', with_Sigma=True, with_dc=True, broadening=None,
-        if self.use_rotations:
+                      transform_to_solver_blocks=True, show_warnings=True):
-            for icrsh in range(self.n_corr_shells):
-                for bname, gf in SK_Sigma_imp[icrsh]:
-                    gf << self.rotloc(icrsh, gf, direction='toGlobal')
-
-    def extract_G_loc(self, mu=None, iw_or_w='iw', with_Sigma=True, with_dc=True, broadening=None):
         r"""
         Extracts the local downfolded Green function by the Brillouin-zone integration of the lattice Green's function.
 
@@ -643,13 +679,20 @@ class SumkDFT(object):
             Imaginary shift for the axis along which the real-axis GF is calculated.
             If not provided, broadening will be set to double of the distance between mesh points in 'mesh'.
             Only relevant for real-frequency GF.
+        transform_to_solver_blocks : bool, optional
+            If True (default), the returned G_loc will be transformed to the block structure ``gf_struct_solver``,
+            else it will be in ``gf_struct_sumk``.
+        show_warnings : bool, optional
+            Displays warning messages during transformation
+            (Only effective if transform_to_solver_blocks = True 
 
         Returns
         -------
-        G_loc_inequiv : list of BlockGf (Green's function) objects
+        G_loc : list of BlockGf (Green's function) objects
-                        List of the local Green's functions for all inequivalent correlated shells, 
+            List of the local Green's functions for all (inequivalent) correlated shells,
             rotated into the corresponding local frames.
-
+            If ``transform_to_solver_blocks`` is True, it will be one per correlated shell, else one per
+            inequivalent correlated shell.
         """
 
         if mu is None:
@@ -708,20 +751,53 @@ class SumkDFT(object):
                     G_loc[icrsh][bname] << self.rotloc(
                         icrsh, gf, direction='toLocal')
 
+        if transform_to_solver_blocks:
+            return self.transform_to_solver_blocks(G_loc, show_warnings=show_warnings)
+
+        return G_loc
+
+    def transform_to_solver_blocks(self, G_loc, G_out=None, show_warnings = True):
+        """ transform G_loc from sumk to solver space
+
+        Parameters
+        ----------
+        G_loc : list of BlockGf
+            a list of one BlockGf per correlated shell with a structure
+            according to ``gf_struct_sumk``, e.g. as returned by
+            :py:meth:`.extract_G_loc` with ``transform_to_solver_blocks=False``.
+        G_out : list of BlockGf
+            a list of one BlockGf per *inequivalent* correlated shell
+            with a structure according to ``gf_struct_solver``.
+            The output Green's function (if not given, a new one is
+            created)
+
+        Returns
+        -------
+        G_out
+        """
+
+        assert isinstance(G_loc, list), "G_loc must be a list (with elements for each correlated shell)"
+
+        if G_out is None:
+            G_out = [self.block_structure.create_gf(ish=ish, mesh=G_loc[self.inequiv_to_corr[ish]].mesh)
+                     for ish in range(self.n_inequiv_shells)]
+        else:
+            for ish in range(self.n_inequiv_shells):
+                self.block_structure.check_gf(G_out, ish=ish)
+
         # transform to CTQMC blocks:
         for ish in range(self.n_inequiv_shells):
-            for block, inner in self.gf_struct_solver[ish].iteritems():
+            self.block_structure.convert_gf(
-                for ind1 in inner:
+                G=G_loc[self.inequiv_to_corr[ish]],
-                    for ind2 in inner:
+                G_struct=None,
-                        block_sumk, ind1_sumk = self.solver_to_sumk[
+                ish_from=self.inequiv_to_corr[ish],
-                            ish][(block, ind1)]
+                ish_to=ish,
-                        block_sumk, ind2_sumk = self.solver_to_sumk[
+                space_from='sumk',
-                            ish][(block, ind2)]
+                G_out=G_out[ish],
-                        G_loc_inequiv[ish][block][ind1, ind2] << G_loc[
+                show_warnings = show_warnings)
-                            self.inequiv_to_corr[ish]][block_sumk][ind1_sumk, ind2_sumk]
 
         # return only the inequivalent shells:
-        return G_loc_inequiv
+        return G_out
 
     def analyse_block_structure(self, threshold=0.00001, include_shells=None, dm=None, hloc=None):
         r"""
@@ -864,7 +940,7 @@ class SumkDFT(object):
             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):
+        if all(isinstance(g_sh.mesh, MeshImFreq) 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],
                 make_copies=False) for g_sh in G]
@@ -872,15 +948,15 @@ class SumkDFT(object):
                 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):
+        elif all(isinstance(g_sh.mesh, MeshImTime) 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):
+        elif all(isinstance(g_sh.mesh, MeshReFreq) 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):
+        elif all(isinstance(g_sh.mesh, MeshReTime) for g_sh in G):
             def get_delta_from_mesh(mesh):
                 w0 = None
                 for w in mesh:
@@ -936,6 +1012,8 @@ class SumkDFT(object):
             the Green's function transformed into the new block structure
         """
 
+        assert isinstance(G, list), "G must be a list (with elements for each correlated shell)"
+
         gf = self._get_hermitian_quantity_from_gf(G)
 
         # initialize the variables
@@ -1004,13 +1082,13 @@ class SumkDFT(object):
         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)
+                for ish in range(self.n_inequiv_shells)],self.corr_to_inequiv)
         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()))
+                gf_function=type(G[ish]._first()), space_from='sumk', space_to='solver')
             for ish in range(self.n_inequiv_shells)]
-
+        #print 'c'
         if analyse_deg_shells:
             self.analyse_deg_shells(G_transformed, threshold, include_shells)
         return G_transformed
@@ -1252,6 +1330,82 @@ class SumkDFT(object):
                         # a block was found, break out of the loop
                         break
     
+    def calculate_diagonalization_matrix(self, prop_to_be_diagonal='eal', calc_in_solver_blocks=False, write_to_blockstructure = True, ish=0):
+        """
+        Calculates the diagonalisation matrix, and (optionally) stores it in the BlockStructure.
+
+        Parameters
+        ----------
+        prop_to_be_diagonal : string, optional
+                              Defines the property to be diagonalized.
+
+                              - 'eal' : local hamiltonian (i.e. crystal field)
+                              - 'dm' : local density matrix
+
+        calc_in_solver_blocks : bool, optional
+                                Whether the property shall be diagonalized in the
+                                full sumk structure, or just in the solver structure.
+                                
+        write_to_blockstructure : bool, optional
+                                Whether the diagonalization matrix shall be written to
+                                the BlockStructure directly.
+        ish : int, optional
+              Number of the correlated shell to be diagonalized.
+
+        Returns
+        -------
+        trafo : dict
+               The transformation matrix for each spin-block in the correlated shell
+
+        """
+        trafo = {}
+        
+        
+        if prop_to_be_diagonal == 'eal':
+            prop = self.eff_atomic_levels()[ish]
+        elif prop_to_be_diagonal == 'dm':
+            prop = self.density_matrix(method='using_point_integration')[ish]
+        else:
+            mpi.report(
+                "calculate_diagonalization_matrix: not a valid quantitiy to be diagonal. Choices are 'eal' or 'dm'.")
+            return 0
+
+        if calc_in_solver_blocks:
+            trafo_tmp = self.block_structure.transformation
+            self.block_structure.transformation = None
+
+            prop_solver = self.block_structure.convert_matrix(prop, space_from='sumk', space_to='solver')
+            t= {}
+            for name in prop_solver:
+                t[name] = numpy.linalg.eigh(prop_solver[name])[1].conjugate().transpose()
+            trafo = self.block_structure.convert_matrix(t, space_from='solver', space_to='sumk')
+            #self.T = numpy.dot(self.T.transpose().conjugate(),
+             #                  self.w).conjugate().transpose()
+            self.block_structure.transformation = trafo_tmp
+        else:
+            for name in prop:
+                t = numpy.linalg.eigh(prop[name])[1].conjugate().transpose()
+                trafo[name] = t
+                # calculate new Transformation matrix
+                #self.T = numpy.dot(self.T.transpose().conjugate(),
+                #                   self.w).conjugate().transpose()
+
+        # measure for the 'unity' of the transformation:
+        #wsqr = sum(abs(self.w.diagonal())**2) / self.w.diagonal().size
+        #return wsqr
+        
+        if write_to_blockstructure:
+            if self.block_structure.transformation == None:
+                self.block_structure.transformation = [{} for icrsh in range(self.n_corr_shells)]
+                for icrsh in range(self. n_corr_shells):
+                    for sp in self.spin_block_names[self.corr_shells[icrsh]['SO']]:
+                        self.block_structure.transformation[icrsh][sp] = numpy.eye(self.corr_shells[icrsh]['dim'], dtype=numpy.complex_)
+
+        
+            self.block_structure.transformation[ish] = trafo
+            
+        return trafo
+        
 
     def density_matrix(self, method='using_gf', beta=40.0):
         """Calculate density matrices in one of two ways.
@@ -1457,9 +1611,10 @@ class SumkDFT(object):
         self.dc_imp = dc_imp
         self.dc_energ = dc_energ
 
-    def calc_dc(self, dens_mat, orb=0, U_interact=None, J_hund=None, use_dc_formula=0, use_dc_value=None):
+    def calc_dc(self, dens_mat, orb=0, U_interact=None, J_hund=None,
+                use_dc_formula=0, use_dc_value=None, transform=True):
         r"""
-        Calculates and sets the double counting corrections.
+        Calculate and set the double counting corrections.
 
         If 'use_dc_value' is provided the double-counting term is uniformly initialized
         with this constant and 'U_interact' and 'J_hund' are ignored.
@@ -1494,7 +1649,9 @@ class SumkDFT(object):
         use_dc_value : float, optional
             Value of the double-counting correction. If specified
             `U_interact`, `J_hund` and `use_dc_formula` are ignored.
-
+        transform : bool
+            whether or not to use the transformation in block_structure
+            to transform the dc
         """
 
         for icrsh in range(self.n_corr_shells):
@@ -1575,6 +1732,11 @@ class SumkDFT(object):
                 mpi.report(
                     "DC for shell %(icrsh)i = %(use_dc_value)f" % locals())
                 mpi.report("DC energy = %s" % self.dc_energ[icrsh])
+            if transform:
+                for sp in spn:
+                    T = self.block_structure.effective_transformation_sumk[icrsh][sp]
+                    self.dc_imp[icrsh][sp] = numpy.dot(T.conjugate().transpose(),
+                            numpy.dot(self.dc_imp[icrsh][sp], T))
 
     def add_dc(self, iw_or_w="iw"):
         r"""
@@ -1606,7 +1768,7 @@ class SumkDFT(object):
 
         return sigma_minus_dc
 
-    def symm_deg_gf(self, gf_to_symm, orb):
+    def symm_deg_gf(self, gf_to_symm, ish=0):
         r"""
         Averages a GF over degenerate shells.
 
@@ -1618,8 +1780,8 @@ class SumkDFT(object):
         ----------
         gf_to_symm : gf_struct_solver like
                      Input and output GF (i.e., it gets overwritten)
-        orb : int
+        ish : int
-              Index of an inequivalent shell.
+              Index of an inequivalent shell. (default value 0)
 
         """
 
@@ -1627,7 +1789,7 @@ class SumkDFT(object):
         # an h5 file, self.deg_shells might be None
         if self.deg_shells is None: return
 
-        for degsh in self.deg_shells[orb]:
+        for degsh in self.deg_shells[ish]:
             # ss will hold the averaged orbitals in the basis where the
             # blocks are all equal
             # i.e. maybe_conjugate(v^dagger gf v)
@@ -1722,8 +1884,10 @@ class SumkDFT(object):
         # collect data from mpi:
         dens = mpi.all_reduce(mpi.world, dens, lambda x, y: x + y)
         mpi.barrier()
-
+        import numpy as np
-        return dens
+        if np.abs(np.imag(dens)) > 1e-20:
+            mpi.report("Warning: Imaginary part in density will be ignored ({})".format(str(np.abs(np.imag(dens)))))
+        return np.real(dens)
 
     def set_mu(self, mu):
         r"""
@@ -2052,3 +2216,31 @@ class SumkDFT(object):
     def __set_deg_shells(self,value):
         self.block_structure.deg_shells = value
     deg_shells = property(__get_deg_shells,__set_deg_shells)
+
+    @property
+    def gf_struct_solver_list(self):
+        return self.block_structure.gf_struct_solver_list
+
+    @property
+    def gf_struct_sumk_list(self):
+        return self.block_structure.gf_struct_sumk_list
+
+    @property
+    def gf_struct_solver_dict(self):
+        return self.block_structure.gf_struct_solver_dict
+
+    @property
+    def gf_struct_sumk_dict(self):
+        return self.block_structure.gf_struct_sumk_dict
+
+    def __get_corr_to_inequiv(self):
+        return self.block_structure.corr_to_inequiv
+    def __set_corr_to_inequiv(self, value):
+        self.block_structure.corr_to_inequiv = value
+    corr_to_inequiv = property(__get_corr_to_inequiv, __set_corr_to_inequiv)
+
+    def __get_inequiv_to_corr(self):
+        return self.block_structure.inequiv_to_corr
+    def __set_inequiv_to_corr(self, value):
+        self.block_structure.inequiv_to_corr = value
+    inequiv_to_corr = property(__get_inequiv_to_corr, __set_inequiv_to_corr)

python/trans_basis.py

@@ -48,7 +48,7 @@ class TransBasis:
         self.T = copy.deepcopy(self.SK.T[0])
         self.w = numpy.identity(SK.corr_shells[0]['dim'])
 
-    def calculate_diagonalisation_matrix(self, prop_to_be_diagonal='eal'):
+    def calculate_diagonalisation_matrix(self, prop_to_be_diagonal='eal', calc_in_solver_blocks = False):
         """
         Calculates the diagonalisation matrix w, and stores it as member of the class.
 
@@ -60,6 +60,10 @@ class TransBasis:
                               - 'eal' : local hamiltonian (i.e. crystal field)
                               - 'dm' : local density matrix
 
+        calc_in_solver_blocks : bool, optional
+                                Whether the property shall be diagonalized in the
+                                full sumk structure, or just in the solver structure.
+
         Returns
         -------
         wsqr : double
@@ -76,6 +80,19 @@ class TransBasis:
                 "trans_basis: not a valid quantitiy to be diagonal. Choices are 'eal' or 'dm'.")
             return 0
 
+        if calc_in_solver_blocks:
+            trafo = self.SK.block_structure.transformation
+            self.SK.block_structure.transformation = None
+
+            prop_solver = self.SK.block_structure.convert_matrix(prop, space_from='sumk', space_to='solver')
+            v= {}
+            for name in prop_solver:
+                v[name] = numpy.linalg.eigh(prop_solver[name])[1]
+            self.w = self.SK.block_structure.convert_matrix(v, space_from='solver', space_to='sumk')['ud' if self.SK.SO else 'up']
+            self.T = numpy.dot(self.T.transpose().conjugate(),
+                               self.w).conjugate().transpose()
+            self.SK.block_structure.transformation = trafo
+        else:
             if self.SK.SO == 0:
                 self.eig, self.w = numpy.linalg.eigh(prop['up'])
                 # calculate new Transformation matrix

test/analyse_block_structure_from_gf.py

@@ -26,7 +26,6 @@ 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
@@ -163,9 +162,9 @@ for conjugate in conjugate_values:
     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_noisy = SK.block_structure.convert_gf(G_noisy, block_structure1, beta = G_noisy.mesh.beta, space_from='sumk')
-    G_pre_transform = SK.block_structure.convert_gf(G_pre_transform, 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, space_from='sumk')
-    G_noisy_pre_transform = SK.block_structure.convert_gf(G_noisy_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, space_from='sumk')
 
     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"

test/blockstructure.py

@@ -2,67 +2,219 @@ from triqs_dft_tools.sumk_dft import *
 from pytriqs.utility.h5diff import h5diff
 from pytriqs.gf import *
 from pytriqs.utility.comparison_tests import assert_block_gfs_are_close
+from scipy.linalg import expm
 from triqs_dft_tools.block_structure import BlockStructure
+import numpy as np
+from pytriqs.utility.h5diff import compare, failures
 
-SK = SumkDFT('blockstructure.in.h5',use_dft_blocks=True)
+
+def cmp(a, b, precision=1.e-15):
+    compare('', a, b, 0, precision)
+    if failures:
+        raise AssertionError('\n'.join(failures))
+
+SK = SumkDFT('blockstructure.in.h5', use_dft_blocks=True)
 
 original_bs = SK.block_structure
+cmp(original_bs.effective_transformation_sumk,
+    [{'down': np.array([[1., 0., 0.],
+                        [0., 1., 0.],
+                        [0., 0., 1.]]),
+      'up': np.array([[1., 0., 0.],
+                      [0., 1., 0.],
+                      [0., 0., 1.]])}])
+cmp(original_bs.effective_transformation_solver,
+    [{'up_0': np.array([[1., 0., 0.],
+                        [0., 1., 0.]]),
+      'up_1': np.array([[0., 0., 1.]]),
+      'down_1': np.array([[0., 0., 1.]]),
+      'down_0': np.array([[1., 0., 0.],
+                          [0., 1., 0.]])}])
+
+created_matrix = original_bs.create_matrix()
+cmp(created_matrix,
+    {'up_0': np.array([[0. + 0.j, 0. + 0.j],
+                       [0. + 0.j, 0. + 0.j]]),
+     'up_1': np.array([[0. + 0.j]]),
+     'down_1': np.array([[0. + 0.j]]),
+     'down_0': np.array([[0. + 0.j, 0. + 0.j],
+                         [0. + 0.j, 0. + 0.j]])})
+
 
 # check pick_gf_struct_solver
 pick1 = original_bs.copy()
 pick1.pick_gf_struct_solver([{'up_0': [1], 'up_1': [0], 'down_1': [0]}])
 
+cmp(pick1.effective_transformation_sumk,
+    [{'down': np.array([[0., 0., 0.],
+                        [0., 0., 0.],
+                        [0., 0., 1.]]),
+      'up': np.array([[0., 0., 0.],
+                      [0., 1., 0.],
+                      [0., 0., 1.]])}])
+cmp(pick1.effective_transformation_solver,
+    [{'up_0': np.array([[0., 1., 0.]]),
+      'up_1': np.array([[0., 0., 1.]]),
+      'down_1': np.array([[0., 0., 1.]])}])
+
 # check loading a block_structure from file
-SK.block_structure = SK.load(['block_structure'],'mod')[0]
+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,
+sk_pick1 = BlockStructure(gf_struct_sumk=SK.gf_struct_sumk,
-                          gf_struct_solver = SK.gf_struct_solver,
+                          gf_struct_solver=SK.gf_struct_solver,
-                          solver_to_sumk = SK.solver_to_sumk,
+                          solver_to_sumk=SK.solver_to_sumk,
-                          sumk_to_solver = SK.sumk_to_solver,
+                          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)
+                          deg_shells=SK.deg_shells,
+                          corr_to_inequiv=SK.corr_to_inequiv)
 assert sk_pick1 == pick1, 'constructing block structure from SumkDFT properties failed'
 
+cmp(pick1.effective_transformation_sumk,
+    [{'down': np.array([[0., 0., 0.],
+                        [0., 0., 0.],
+                        [0., 0., 1.]]),
+      'up': np.array([[0., 0., 0.],
+                      [0., 1., 0.],
+                      [0., 0., 1.]])}])
+cmp(pick1.effective_transformation_solver,
+    [{'up_0': np.array([[0., 1., 0.]]),
+      'up_1': np.array([[0., 0., 1.]]),
+      'down_1': np.array([[0., 0., 1.]])}])
+
 # check pick_gf_struct_sumk
 pick2 = original_bs.copy()
-pick2.pick_gf_struct_sumk([{'up': [1, 2], 'down': [0,1]}])
+pick2.pick_gf_struct_sumk([{'up': [1, 2], 'down': [0, 1]}])
+
+cmp(pick2.effective_transformation_sumk,
+    [{'down': np.array([[1., 0., 0.],
+                        [0., 1., 0.],
+                        [0., 0., 0.]]),
+      'up': np.array([[0., 0., 0.],
+                      [0., 1., 0.],
+                      [0., 0., 1.]])}])
+cmp(pick2.effective_transformation_solver,
+    [{'up_0': np.array([[0., 1., 0.]]),
+      'up_1': np.array([[0., 0., 1.]]),
+      'down_0': np.array([[1., 0., 0.],
+                          [0., 1., 0.]])}])
+
+pick3 = pick2.copy()
+pick3.transformation = [np.reshape(range(9), (3, 3))]
+cmp(pick3.effective_transformation_sumk,
+    [{'down': np.array([[0, 1, 2],
+                        [3, 4, 5],
+                        [0, 0, 0]]),
+      'up': np.array([[0, 0, 0],
+                      [3, 4, 5],
+                      [6, 7, 8]])}])
+cmp(pick3.effective_transformation_solver,
+    [{'up_0': np.array([[3, 4, 5]]),
+      'up_1': np.array([[6, 7, 8]]),
+      'down_0': np.array([[0, 1, 2],
+                          [3, 4, 5]])}])
+
+pick4 = original_bs.copy()
+pick4.transformation = [np.array([[0, 1, 0], [1, 0, 0], [0, 0, 1]])]
+pick4.pick_gf_struct_sumk([{'up': [1, 2], 'down': [0, 1]}])
+cmp(pick2.gf_struct_sumk, pick4.gf_struct_sumk)
+cmp(pick2.gf_struct_solver, pick4.gf_struct_solver)
+assert pick4.sumk_to_solver == [{('up', 0): ('up_0', 0),
+                                 ('up', 1): (None, None),
+                                 ('up', 2): ('up_1', 0),
+                                 ('down', 2): (None, None),
+                                 ('down', 1): ('down_0', 1),
+                                 ('down', 0): ('down_0', 0)}]
+assert pick4.solver_to_sumk == [{('up_1', 0): ('up', 2),
+                                 ('up_0', 0): ('up', 0),
+                                 ('down_0', 0): ('down', 0),
+                                 ('down_0', 1): ('down', 1)}]
 
 # check map_gf_struct_solver
-mapping = [{ ('down_0', 0):('down', 0),
+mapping = [{('down_0', 0): ('down', 0),
-             ('down_0', 1):('down', 2),
+            ('down_0', 1): ('down', 2),
-             ('down_1', 0):('down', 1),
+            ('down_1', 0): ('down', 1),
-             ('up_0', 0)  :('down_1', 0),
+            ('up_0', 0): ('down_1', 0),
-             ('up_0', 1)  :('up_0', 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)
+G1 = original_bs.create_gf(beta=40, n_points=3)
-i = 1
+widths = dict(up_0=1, up_1=2, down_0=4, down_1=3)
-for block,gf in G1:
+for block, gf in G1:
-    gf << SemiCircular(i)
+    gf << SemiCircular(widths[block])
-    i+=1
+original_bs.check_gf(G1)
+original_bs.check_gf([G1])
 
 # 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)
+import warnings
+with warnings.catch_warnings(record=True) as w:
+    G2 = map1.convert_gf(G1, original_bs, beta=40, n_points=3,
+                         show_warnings=True)
+    assert len(w) == 1
+    assert issubclass(w[-1].category, UserWarning)
+    assert "Block up_1 maximum difference" in str(w[-1].message)
+
+m2 = map1.convert_matrix(created_matrix, original_bs, show_warnings=True)
+cmp(m2,
+    {'down': np.array([[0. + 0.j, 0. + 0.j, 0. + 0.j],
+                       [0. + 0.j, 0. + 0.j, 0. + 0.j],
+                       [0. + 0.j, 0. + 0.j, 0. + 0.j]]),
+     'up_0': np.array([[0. + 0.j]]),
+     'down_1': np.array([[0. + 0.j]])})
 
 # check full_structure
-full = BlockStructure.full_structure([{'up_0': [0, 1], 'up_1': [0], 'down_1': [0], 'down_0': [0, 1]}],None)
+full = BlockStructure.full_structure(
+    [{'up_0': [0, 1], 'up_1': [0], 'down_1': [0], 'down_0': [0, 1]}], None)
+
+G_sumk = BlockGf(mesh=G1.mesh, gf_struct=original_bs.gf_struct_sumk[0])
+for i in range(3):
+    G_sumk['up'][i, i] << SemiCircular(1 if i < 2 else 2)
+    G_sumk['down'][i, i] << SemiCircular(4 if i < 2 else 3)
+G3 = original_bs.convert_gf(G_sumk,
+                            None,
+                            space_from='sumk',
+                            beta=40,
+                            n_points=3)
+assert_block_gfs_are_close(G1, G3)
+
+# check convert_gf with transformation
+# np.random.seed(894892309)
+H = np.random.rand(3, 3) + 1.0j * np.random.rand(3, 3)
+H = H + H.conjugate().transpose()
+T = expm(1.0j * H)
+G_T = G_sumk.copy()
+for block, gf in G_T:
+    gf.from_L_G_R(T.conjugate().transpose(), gf, T)
+transformed_bs = original_bs.copy()
+transformed_bs.transformation = [T]
+G_bT = transformed_bs.convert_gf(G_T, None, space_from='sumk',
+                                 beta=40, n_points=3)
+assert_block_gfs_are_close(G1, G_bT)
+
+assert original_bs.gf_struct_sumk_list ==\
+    [[('up', [0, 1, 2]), ('down', [0, 1, 2])]]
+assert original_bs.gf_struct_solver_dict ==\
+    [{'up_0': [0, 1], 'up_1': [0], 'down_1': [0], 'down_0': [0, 1]}]
+assert original_bs.gf_struct_sumk_dict ==\
+    [{'down': [0, 1, 2], 'up': [0, 1, 2]}]
+assert original_bs.gf_struct_solver_list ==\
+    [[('down_0', [0, 1]), ('down_1', [0]), ('up_0', [0, 1]), ('up_1', [0])]]
 
 # check __eq__
-assert full==full, 'equality not correct (equal structures not equal)'
+assert full == full, 'equality not correct (equal structures not equal)'
-assert pick1==pick1, '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 pick1 != pick2, 'equality not correct (different structures not different)'
-assert original_bs!=offd, '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:
+    with HDFArchive('blockstructure.out.h5', 'w') as ar:
         ar['original_bs'] = original_bs
         ar['pick1'] = pick1
         ar['pick2'] = pick2
@@ -75,10 +227,10 @@ if mpi.is_master_node():
     # 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,\
+    with HDFArchive('blockstructure.out.h5', 'r') as ar,\
-         HDFArchive('blockstructure.ref.h5','r') as ar2:
+            HDFArchive('blockstructure.ref.h5', 'r') as ar2:
         for k in ar2:
-                if isinstance(ar[k],BlockGf):
+            if isinstance(ar[k], BlockGf):
-                    assert_block_gfs_are_close(ar[k],ar2[k],1.e-6)
+                assert_block_gfs_are_close(ar[k], ar2[k], 1.e-6)
             else:
-                    assert ar[k]==ar2[k], '{} not equal'.format(k)
+                assert ar[k] == ar2[k], '{} not equal'.format(k)