SumkDFT: analyze_block_structure_from_gf for Gf Re/Im Time/Freq

python/sumk_dft.py

@@ -883,7 +883,7 @@ class SumkDFT(object):
             the Green's function transformed into the new block structure
         """
 
-        # make a GfImTime from the supplied G
+        # make a GfImTime from the supplied GfImFreq
         if all(isinstance(g_sh._first(), GfImFreq) for g_sh in G):
             gf = [BlockGf(name_block_generator = [(name, GfImTime(beta=block.mesh.beta,
                 indices=block.indices,n_points=len(block.mesh)+1)) for name, block in g_sh])
@@ -891,9 +891,37 @@ class SumkDFT(object):
             for ish in range(len(gf)):
                 for name, g in gf[ish]:
                     g.set_from_inverse_fourier(G[ish][name])
-        else:
+        # keep a GfImTime from the supplied GfImTime
-            assert all(isinstance(g_sh._first(), GfImTime) for g_sh in G), "G must be a BlockGf of either GfImFreq or GfImTime"
+        elif all(isinstance(g_sh._first(), GfImTime) for g_sh in G):
             gf = G
+        # make a spectral function from the supplied GfReFreq
+        elif all(isinstance(g_sh._first(), GfReFreq) for g_sh in G):
+            gf = [g_sh.copy() for g_sh in G]
+            for ish in range(len(gf)):
+                for name, g in gf[ish]:
+                    g << 1.0j*(g-g.conjugate().transpose())/2.0/numpy.pi
+        elif all(isinstance(g_sh._first(), GfReTime) for g_sh in G):
+            def get_delta_from_mesh(mesh):
+                w0 = None
+                for w in mesh:
+                    if w0 is None:
+                        w0 = w
+                    else:
+                        return w-w0
+            gf = [BlockGf(
+                name_block_generator = [(name,
+                    GfReFreq(
+                        window=(-numpy.pi*(len(block.mesh)-1) / (len(block.mesh)*get_delta_from_mesh(block.mesh)), numpy.pi*(len(block.mesh)-1) / (len(block.mesh)*get_delta_from_mesh(block.mesh))),
+                        n_points=len(block.mesh),
+                        indices=block.indices)) for name, block in g_sh])
+                for g_sh in G]
+
+            for ish in range(len(gf)):
+                for name, g in gf[ish]:
+                    g.set_from_fourier(G[ish][name])
+                    g << 1.0j*(g-g.conjugate().transpose())/2.0/numpy.pi
+        else:
+            raise Exception("G must be a list of BlockGf of either GfImFreq, GfImTime, GfReFreq or GfReTime")
 
         # initialize the variables
         self.gf_struct_solver = [{} for ish in range(self.n_inequiv_shells)]
@@ -964,7 +992,8 @@ class SumkDFT(object):
                 for ish in range(self.n_inequiv_shells)],None)
         G_transformed = [
             self.block_structure.convert_gf(G[ish],
-                full_structure, ish, beta=G[ish].mesh.beta, show_warnings=threshold)
+                full_structure, ish, mesh=G[ish].mesh.copy(), show_warnings=threshold,
+                gf_function=type(G[ish]._first()))
             for ish in range(self.n_inequiv_shells)]
 
         if analyse_deg_shells:
@@ -1002,7 +1031,7 @@ class SumkDFT(object):
             null_space = compress(null_mask, vh, axis=0)
             return null_space.conjugate().transpose()
 
-        # make a GfImTime from the supplied G
+        # make a GfImTime from the supplied GfImFreq
         if all(isinstance(g_sh._first(), GfImFreq) for g_sh in G):
             gf = [BlockGf(name_block_generator = [(name, GfImTime(beta=block.mesh.beta,
                 indices=block.indices,n_points=len(block.mesh)+1)) for name, block in g_sh])
@@ -1010,9 +1039,37 @@ class SumkDFT(object):
             for ish in range(len(gf)):
                 for name, g in gf[ish]:
                     g.set_from_inverse_fourier(G[ish][name])
-        else:
+        # keep a GfImTime from the supplied GfImTime
-            assert all(isinstance(g_sh._first(), GfImTime) for g_sh in G), "G must be a BlockGf of either GfImFreq or GfImTime"
+        elif all(isinstance(g_sh._first(), GfImTime) for g_sh in G):
             gf = G
+        # make a spectral function from the supplied GfReFreq
+        elif all(isinstance(g_sh._first(), GfReFreq) for g_sh in G):
+            gf = [g_sh.copy() for g_sh in G]
+            for ish in range(len(gf)):
+                for name, g in gf[ish]:
+                    g << 1.0j*(g-g.conjugate().transpose())/2.0/numpy.pi
+        elif all(isinstance(g_sh._first(), GfReTime) for g_sh in G):
+            def get_delta_from_mesh(mesh):
+                w0 = None
+                for w in mesh:
+                    if w0 is None:
+                        w0 = w
+                    else:
+                        return w-w0
+            gf = [BlockGf(
+                name_block_generator = [(name,
+                    GfReFreq(
+                        window=(-numpy.pi*(len(block.mesh)-1) / (len(block.mesh)*get_delta_from_mesh(block.mesh)), numpy.pi*(len(block.mesh)-1) / (len(block.mesh)*get_delta_from_mesh(block.mesh))),
+                        n_points=len(block.mesh),
+                        indices=block.indices)) for name, block in g_sh])
+                for g_sh in G]
+
+            for ish in range(len(gf)):
+                for name, g in gf[ish]:
+                    g.set_from_fourier(G[ish][name])
+                    g << 1.0j*(g-g.conjugate().transpose())/2.0/numpy.pi
+        else:
+            raise Exception("G must be a list of BlockGf of either GfImFreq, GfImTime, GfReFreq or GfReTime")
 
         if include_shells is None:
             # include all shells
@@ -1606,7 +1663,7 @@ class SumkDFT(object):
                 # the helper is in the basis where the blocks are all equal
                 helper.from_L_G_R(v.conjugate().transpose(), gf_to_symm[key], v)
                 if C:
-                    conjugate_in_tau(helper, in_place=True)
+                    helper << helper.transpose()
                 # average over all shells
                 ss += helper / (1.0 * n_deg)
             # now put back the averaged gf to all shells
@@ -1618,11 +1675,7 @@ class SumkDFT(object):
                     v = numpy.eye(*ss.target_shape)
                     C = False
                 if C:
-                    # we could also use
+                    gf_to_symm[key].from_L_G_R(v, ss.transpose(), v.conjugate().transpose())
-                    # gf_to_symm[key].from_L_G_R(v, conjugate_in_tau(ss), v.conjugate().transpose())
-                    # but this is less memory-intensive:
-                    gf_to_symm[key].from_L_G_R(v.conjugate(), ss, v.transpose())
-                    conjugate_in_tau(gf_to_symm[key], in_place=True)
                 else:
                     gf_to_symm[key].from_L_G_R(v, ss, v.conjugate().transpose())
 
@@ -2015,38 +2068,3 @@ class SumkDFT(object):
     def __set_deg_shells(self,value):
         self.block_structure.deg_shells = value
     deg_shells = property(__get_deg_shells,__set_deg_shells)
-
-# a helper function
-def conjugate_in_tau(gf_im_freq, in_place=False):
-    """ Calculate the conjugate in tau of a GfImFreq
-
-    Parameters
-    ----------
-    gf_im_freq : GfImFreq of BlockGf
-        the Green's function
-    in_place : whether to modify the gf_im_freq object (True) or return a copy (False)
-
-    Returns
-    -------
-    ret : GfImFreq of BlockGf
-        the Green's function that has been FT to G(tau), conjugated, and
-        FT back
-    """
-    if in_place:
-        ret = gf_im_freq
-    else:
-        ret = gf_im_freq.copy()
-    if isinstance(ret, BlockGf):
-        for name, gf in ret:
-            conjugate_in_tau(gf, in_place=True)
-    else:
-        """ there is an easier way to do this, namely to make
-            ret.data[:,:,:] = gf_im_freq.data[::-1,:,:].conjugate()
-            ret.tail.data[:,:,:] = gf_im_freq.tail.data.conjugate()
-            but this relies on symmetric Matsubara meshes and is maybe
-            not safe enough"""
-        G_tau = GfImTime(beta=gf_im_freq.mesh.beta,
-                    indices=gf_im_freq.indices,n_points=len(gf_im_freq.mesh)+1)
-        G_tau.set_from_inverse_fourier(gf_im_freq)
-        ret.set_from_fourier(G_tau.conjugate())
-    return ret

test/CMakeLists.txt

@@ -5,7 +5,7 @@ file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/${all_h5_files} DESTINATION ${CMAKE_CURREN
 FILE(COPY SrVO3.pmat SrVO3.struct SrVO3.outputs SrVO3.oubwin SrVO3.ctqmcout SrVO3.symqmc SrVO3.sympar SrVO3.parproj SrIrO3_rot.h5 hk_convert_hamiltonian.hk LaVO3-Pnma_hr.dat LaVO3-Pnma.inp DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
 
 # List all tests
-set(all_tests wien2k_convert hk_convert w90_convert sumkdft_basic srvo3_Gloc srvo3_transp sigma_from_file blockstructure analyze_block_structure_from_gf)
+set(all_tests wien2k_convert hk_convert w90_convert sumkdft_basic srvo3_Gloc srvo3_transp sigma_from_file blockstructure analyze_block_structure_from_gf analyze_block_structure_from_gf2)
 
 foreach(t ${all_tests})
   add_test(NAME ${t} COMMAND python ${CMAKE_CURRENT_SOURCE_DIR}/${t}.py)

test/analyze_block_structure_from_gf.py

@@ -1,5 +1,5 @@
 from pytriqs.gf import *
-from sumk_dft import SumkDFT, conjugate_in_tau
+from sumk_dft import SumkDFT
 from scipy.linalg import expm
 import numpy as np
 from pytriqs.utility.comparison_tests import assert_gfs_are_close, assert_arrays_are_close, assert_block_gfs_are_close
@@ -68,7 +68,7 @@ for d in SK.deg_shells[0]:
         normalized_gf = G_new[0][key].copy()
         normalized_gf.from_L_G_R(d[key][0].conjugate().transpose(), G_new[0][key], d[key][0])
         if d[key][1]:
-            conjugate_in_tau(normalized_gf, in_place=True)
+            normalized_gf << normalized_gf.transpose()
         normalized_gfs.append(normalized_gf)
     for i in range(len(normalized_gfs)):
         for j in range(i+1,len(normalized_gfs)):
@@ -162,8 +162,8 @@ for conjugate in conjugate_values:
         G_noisy['ud'][i:i+2,i:i+2].from_L_G_R(T, G_noisy['ud'][i:i+2,i:i+2], T.conjugate().transpose())
         # if that block shall be conjugated, go ahead and do it
         if conjugate[i//2]:
-            conjugate_in_tau(G[0]['ud'][i:i+2,i:i+2], in_place=True)
+            G[0]['ud'][i:i+2,i:i+2] << G[0]['ud'][i:i+2,i:i+2].transpose()
-            conjugate_in_tau(G_noisy['ud'][i:i+2,i:i+2], in_place=True)
+            G_noisy['ud'][i:i+2,i:i+2] << G_noisy['ud'][i:i+2,i:i+2].transpose()
 
     # analyse the block structure
     G_new = SK.analyse_block_structure_from_gf(G)
@@ -196,7 +196,7 @@ for conjugate in conjugate_values:
             normalized_gf = G_new[0][key].copy()
             normalized_gf.from_L_G_R(d[key][0].conjugate().transpose(), G_new[0][key], d[key][0])
             if d[key][1]:
-                conjugate_in_tau(normalized_gf, in_place=True)
+                normalized_gf << normalized_gf.transpose()
             normalized_gfs.append(normalized_gf)
         for i in range(len(normalized_gfs)):
             for j in range(i+1,len(normalized_gfs)):

test/analyze_block_structure_from_gf2.py

@@ -0,0 +1,115 @@
+from pytriqs.gf import *
+from sumk_dft import SumkDFT
+import numpy as np
+from pytriqs.utility.comparison_tests import assert_block_gfs_are_close
+
+# here we test the SK.analyze_block_structure_from_gf function
+# with GfReFreq, GfReTime
+
+
+# helper function to get random Hermitian matrix
+def get_random_hermitian(dim):
+    herm = np.random.rand(dim,dim)+1.0j*np.random.rand(dim,dim)
+    herm = herm + herm.conjugate().transpose()
+    return herm
+
+# helper function to get random unitary matrix
+def get_random_transformation(dim):
+    herm = get_random_hermitian(dim)
+    T = expm(1.0j*herm)
+    return T
+
+# construct a random block-diagonal Hloc
+Hloc = np.zeros((10,10), dtype=np.complex_)
+# the Hloc of the first three 2x2 blocks is equal
+Hloc0 = get_random_hermitian(2)
+Hloc[:2,:2] = Hloc0
+Hloc[2:4,2:4] = Hloc0
+Hloc[4:6,4:6] = Hloc0
+# the Hloc of the last two 2x2 blocks is equal
+Hloc1 = get_random_hermitian(2)
+Hloc[6:8,6:8] = Hloc1
+Hloc[8:,8:] = Hloc1
+# construct the hybridization delta
+# this is equal for all 2x2 blocks
+V = get_random_hermitian(2) # the hopping elements from impurity to bath
+b1 = np.random.rand() # the bath energy of the first bath level
+b2 = np.random.rand() # the bath energy of the second bath level
+delta = GfReFreq(window=(-5,5), indices=range(2), n_points=1001)
+delta[0,0] << (V[0,0]*V[0,0].conjugate()*inverse(Omega-b1)+V[0,1]*V[0,1].conjugate()*inverse(Omega-b2+0.02j))/2.0
+delta[0,1] << (V[0,0]*V[1,0].conjugate()*inverse(Omega-b1)+V[0,1]*V[1,1].conjugate()*inverse(Omega-b2+0.02j))/2.0
+delta[1,0] << (V[1,0]*V[0,0].conjugate()*inverse(Omega-b1)+V[1,1]*V[0,1].conjugate()*inverse(Omega-b2+0.02j))/2.0
+delta[1,1] << (V[1,0]*V[1,0].conjugate()*inverse(Omega-b1)+V[1,1]*V[1,1].conjugate()*inverse(Omega-b2+0.02j))/2.0
+# construct G
+G = BlockGf(name_block_generator=(('ud',GfReFreq(window=(-5,5), indices=range(10), n_points=1001)),))
+for i in range(0,10,2):
+    G['ud'][i:i+2,i:i+2] << inverse(Omega-delta+0.02j)
+G['ud'] << inverse(inverse(G['ud']) - Hloc)
+
+
+SK = SumkDFT(hdf_file = 'SrIrO3_rot.h5', use_dft_blocks=False)
+G_new = SK.analyse_block_structure_from_gf([G])
+G_new_symm = G_new[0].copy()
+SK.symm_deg_gf(G_new_symm, 0)
+assert_block_gfs_are_close(G_new[0], G_new_symm)
+
+
+assert SK.gf_struct_sumk == [[('ud', [0, 1, 2, 3, 4, 5, 6, 7, 8, 9])], [('ud', [0, 1, 2, 3, 4, 5, 6, 7, 8, 9])]],\
+    "wrong gf_struct_sumk"
+for i in range(5):
+    assert 'ud_{}'.format(i) in SK.gf_struct_solver[0], "missing block"
+    assert SK.gf_struct_solver[0]['ud_{}'.format(i)] == range(2), "wrong block size"
+for i in range(10):
+    assert SK.sumk_to_solver[0]['ud',i] == ('ud_{}'.format(i/2), i%2), "wrong mapping"
+
+assert len(SK.deg_shells[0]) == 2, "wrong number of equivalent groups found"
+assert sorted([len(d) for d in SK.deg_shells[0]]) == [2,3], "wrong number of members in the equivalent groups found"
+for d in SK.deg_shells[0]:
+    if len(d)==2:
+        assert 'ud_3' in d, "shell ud_3 missing"
+        assert 'ud_4' in d, "shell ud_4 missing"
+    if len(d)==3:
+        assert 'ud_0' in d, "shell ud_0 missing"
+        assert 'ud_1' in d, "shell ud_1 missing"
+        assert 'ud_2' in d, "shell ud_2 missing"
+
+
+
+def get_delta_from_mesh(mesh):
+    w0 = None
+    for w in mesh:
+        if w0 is None:
+            w0 = w
+        else:
+            return w-w0
+
+Gt = BlockGf(name_block_generator = [(name,
+            GfReTime(window=(-np.pi*(len(block.mesh)-1) / (len(block.mesh)*get_delta_from_mesh(block.mesh)), np.pi*(len(block.mesh)-1) / (len(block.mesh)*get_delta_from_mesh(block.mesh))),
+                n_points=len(block.mesh),
+                indices=block.indices)) for name, block in G])
+
+Gt['ud'].set_from_inverse_fourier(G['ud'])
+
+G_new = SK.analyse_block_structure_from_gf([Gt])
+G_new_symm = G_new[0].copy()
+SK.symm_deg_gf(G_new_symm, 0)
+assert_block_gfs_are_close(G_new[0], G_new_symm)
+
+assert SK.gf_struct_sumk == [[('ud', [0, 1, 2, 3, 4, 5, 6, 7, 8, 9])], [('ud', [0, 1, 2, 3, 4, 5, 6, 7, 8, 9])]],\
+    "wrong gf_struct_sumk"
+for i in range(5):
+    assert 'ud_{}'.format(i) in SK.gf_struct_solver[0], "missing block"
+    assert SK.gf_struct_solver[0]['ud_{}'.format(i)] == range(2), "wrong block size"
+for i in range(10):
+    assert SK.sumk_to_solver[0]['ud',i] == ('ud_{}'.format(i/2), i%2), "wrong mapping"
+
+assert len(SK.deg_shells[0]) == 2, "wrong number of equivalent groups found"
+assert sorted([len(d) for d in SK.deg_shells[0]]) == [2,3], "wrong number of members in the equivalent groups found"
+for d in SK.deg_shells[0]:
+    if len(d)==2:
+        assert 'ud_3' in d, "shell ud_3 missing"
+        assert 'ud_4' in d, "shell ud_4 missing"
+    if len(d)==3:
+        assert 'ud_0' in d, "shell ud_0 missing"
+        assert 'ud_1' in d, "shell ud_1 missing"
+        assert 'ud_2' in d, "shell ud_2 missing"