block_structure: convert_gf with transformation

python/block_structure.py

@@ -29,6 +29,7 @@ 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):
@@ -209,10 +210,23 @@ class BlockStructure(object):
 
     @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.corr_to_inequiv is None:
+            return None
+        ret = []
+        for icrsh, ish in enumerate(self.corr_to_inequiv):
+            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
@@ -671,8 +685,6 @@ class BlockStructure(object):
             options passed to the constructor for the new Gf
         """
 
-        # TODO: use effective_transformation here
-
         if ish is not None:
             warn(
                 'The parameter ish in convert_gf is deprecated. Use ish_from and ish_to instead.')
@@ -691,55 +703,61 @@ class BlockStructure(object):
             G_struct = self
 
         if space_from == 'solver':
-            gf_struct_in = G_struct.gf_struct_solver[ish_from]
+            gf_struct_from = G_struct.gf_struct_solver[ish_from]
+            eff_trans_from = G_struct.effective_transformation_solver[ish_from]
+            block_mapping_from = G_struct.sumk_to_solver_block[ish_from]
         elif space_from == 'sumk':
-            gf_struct_in = G_struct.gf_struct_sumk_dict[ish_from]
+            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'.")
+
         # create a Green's function to hold the result
         if G_out is None:
             G_out = self.create_gf(ish=ish_to, space=space_to, **kwargs)
         else:
             self.check_gf(G_out, ish=ish_to, space=space_to)
 
-        for block in gf_struct_in.keys():
+        for block_to in gf_struct_to.keys():
-            for i1 in gf_struct_in[block]:
+            block_intermediate = block_mapping_to[block_to]
-                for i2 in gf_struct_in[block]:
+            block_from = block_mapping_from[block_intermediate]
-                    if space_from == 'sumk':
+            T_to = eff_trans_to[block_to]
-                        i1_sumk = (block, i1)
+            g_help = G_out[block_to].copy()
-                        i2_sumk = (block, i2)
+            for block in block_from:
-                    elif space_from == 'solver':
+                T_from = eff_trans_from[block]
-                        i1_sumk = G_struct.solver_to_sumk[ish_from][(block, i1)]
+                g_help.from_L_G_R(np.dot(T_to, np.conjugate(np.transpose(T_from))),
-                        i2_sumk = G_struct.solver_to_sumk[ish_from][(block, i2)]
+                                  G[block],
+                                  np.dot(T_from, np.conjugate(np.transpose(T_to))))
+                G_out[block_to] << G_out[block_to] + g_help
 
-                    if space_to == 'sumk':
+        if show_warnings:
-                        i1_sol = i1_sumk
+            # we back-transform it
-                        i2_sol = i2_sumk
+            G_back = G_struct.convert_gf(G_out, self, ish_from=ish_to,
-                    elif space_to == 'solver':
+                                         ish_to=ish_from,
-                        i1_sol = self.sumk_to_solver[ish_to][i1_sumk]
+                                         show_warnings=False,  # else we get an endless loop
-                        i2_sol = self.sumk_to_solver[ish_to][i2_sumk]
+                                         space_from=space_to, space_to=space_from, **kwargs)
-                    else:
+            for name, gf in G_back:
-                        raise Exception(
+                maxdiff = np.max(np.abs(G_back[name].data - G[name].data),
-                            "Argument space_to has to be either 'solver' or 'sumk'.")
+                                 axis=0)
-
+                if np.any(maxdiff > warning_threshold):
-                    if i1_sol[0] is None or i2_sol[0] is None:
+                    warn('Block {} maximum difference:\n'.format(name)
-                        if show_warnings:
+                         + str(maxdiff))
-                            if mpi.is_master_node():
-                                warn(('Element {},{} of block {} of G is not present ' +
-                                      'in the new structure').format(i1, i2, block))
-                        continue
-                    if i1_sol[0] != i2_sol[0]:
-                        if (show_warnings and
-                                np.max(np.abs(G[block][i1, i2].data)) > warning_threshold):
-                            if mpi.is_master_node():
-                                warn(('Element {},{} of block {} of G is approximated ' +
-                                      'to zero to match the new structure. Max abs value: {}').format(
-                                    i1, i2, block, np.max(np.abs(G[block][i1, i2].data))))
-                        continue
-                    G_out[i1_sol[0]][i1_sol[1], i2_sol[1]] = \
-                        G[block][i1, i2]
         return G_out
 
     def approximate_as_diagonal(self):

test/blockstructure.py

@@ -2,6 +2,7 @@ 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
@@ -127,7 +128,14 @@ 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)
+
 
 # check full_structure
 full = BlockStructure.full_structure(
@@ -144,6 +152,20 @@ G3 = original_bs.convert_gf(G_sumk,
                             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 ==\