analyze_block_structure_from_gf test

test/CMakeLists.txt

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

test/analyze_block_structure_from_gf.py

@@ -0,0 +1,179 @@
+from pytriqs.gf import *
+from sumk_dft import SumkDFT, conjugate_in_tau
+from scipy.linalg import expm
+import numpy as np
+from pytriqs.utility.comparison_tests import assert_gfs_are_close, assert_arrays_are_close
+from pytriqs.archive import *
+import itertools
+
+# The full test checks all different possible combinations of conjugated
+# blocks. This takes a few minutes. For a quick test, just checking one
+# random value suffices.
+# (this parameter affects the second test)
+full_test = False
+
+#######################################################################
+# First test                                                          #
+# where we check the analyse_block_structure_from_gf function         #
+# for the SrIrO3_rot.h5 file                                          #
+#######################################################################
+
+beta = 40
+SK = SumkDFT(hdf_file = 'SrIrO3_rot.h5')
+Sigma = SK.block_structure.create_gf(beta=beta)
+SK.put_Sigma([Sigma])
+G = SK.extract_G_loc()
+
+# the original block structure
+block_structure1 = SK.block_structure.copy()
+
+G_new = SK.analyse_block_structure_from_gf(G)
+
+# the new block structure
+block_structure2 = SK.block_structure.copy()
+
+with HDFArchive('analyze_block_structure_from_gf.out.h5','w') as ar:
+    ar['bs1'] = block_structure1
+    ar['bs2'] = block_structure2
+
+# check whether the block structure is the same as in the reference
+with HDFArchive('analyze_block_structure_from_gf.out.h5','r') as ar,\
+     HDFArchive('analyze_block_structure_from_gf.ref.h5','r') as ar2:
+    assert ar['bs1'] == ar2['bs1'], 'bs1 not equal'
+    a1 = ar['bs2']
+    a2 = ar2['bs2']
+    assert a1==block_structure2, "writing/reading block structure incorrect"
+    # we set the deg_shells to None because the transformation matrices
+    # have a phase freedom and will, therefore, not be equal in general
+    a1.deg_shells = None
+    a2.deg_shells = None
+    assert a1==a2, 'bs2 not equal'
+
+# check if deg shells are correct
+assert len(SK.deg_shells[0])==1, "wrong number of equivalent groups"
+
+# check if the Green's functions that are found to be equal in the
+# routine are indeed equal
+for d in SK.deg_shells[0]:
+    assert len(d)==2, "wrong number of shells in equivalent group"
+    # the convention is that for every degenerate shell, the transformation
+    # matrix v and the conjugate bool is saved
+    # then,
+    # maybe_conjugate1( v1^dagger G1 v1 ) = maybe_conjugate2( v2^dagger G2 v2 )
+    # therefore, to test, we calculate
+    #   maybe_conjugate( v^dagger G v )
+    # for all degenerate shells and check that they are all equal
+    normalized_gfs = []
+    for key in d:
+        normalized_gf = G_new[0][key].copy()
+        normalized_gf.from_L_G_R(d[key][0].conjugate().transpose(), G_new[0][key], d[key][0])
+        if d[key][1]:
+            conjugate_in_tau(normalized_gf, in_place=True)
+        normalized_gfs.append(normalized_gf)
+    for i in range(len(normalized_gfs)):
+        for j in range(i+1,len(normalized_gfs)):
+            assert_arrays_are_close(normalized_gfs[i].data, normalized_gfs[j].data, 1.e-5)
+            # the tails have to be compared using a relative error
+            for o in range(normalized_gfs[i].tail.order_min,normalized_gfs[i].tail.order_max+1):
+                if np.abs(normalized_gfs[i].tail[o][0,0]) < 1.e-10:
+                    continue
+                assert np.max(np.abs((normalized_gfs[i].tail[o]-normalized_gfs[j].tail[o])/(normalized_gfs[i].tail[o][0,0]))) < 1.e-5, \
+                        "tails are different"
+
+#######################################################################
+# Second test                                                         #
+# where a Green's function is constructed from a random model         #
+# and the analyse_block_structure_from_gf function is tested for that #
+# model                                                               #
+#######################################################################
+
+# helper function to get random Hermitian matrix
+def get_random_hermitian(dim):
+    herm = np.random.rand(dim,dim)+1.0j*np.random.rand(dim,dim)
+    herm = herm + herm.conjugate().transpose()
+    return herm
+
+# helper function to get random unitary matrix
+def get_random_transformation(dim):
+    herm = get_random_hermitian(dim)
+    T = expm(1.0j*herm)
+    return T
+
+# we will conjugate the Green's function blocks according to the entries
+# of conjugate_values
+# for each of the 5 blocks that will be constructed, there is an entry
+# True or False that says whether it will be conjugated
+if full_test:
+    # in the full test we check all combinations
+    conjugate_values = list(itertools.product([False, True], repeat=5))
+else:
+    # in the quick test we check a random combination
+    conjugate_values = [np.random.rand(5)>0.5]
+
+for conjugate in conjugate_values:
+    # construct a random block-diagonal Hloc
+    Hloc = np.zeros((10,10), dtype=np.complex_)
+    # the Hloc of the first three 2x2 blocks is equal
+    Hloc0 = get_random_hermitian(2)
+    Hloc[:2,:2] = Hloc0
+    Hloc[2:4,2:4] = Hloc0
+    Hloc[4:6,4:6] = Hloc0
+    # the Hloc of the last two 2x2 blocks is equal
+    Hloc1 = get_random_hermitian(2)
+    Hloc[6:8,6:8] = Hloc1
+    Hloc[8:,8:] = Hloc1
+    # construct the hybridization delta
+    # this is equal for all 2x2 blocks
+    V = get_random_hermitian(2) # the hopping elements from impurity to bath
+    b1 = np.random.rand() # the bath energy of the first bath level
+    b2 = np.random.rand() # the bath energy of the second bath level
+    delta = G[0]['ud'][:2,:2].copy()
+    delta[0,0] << (V[0,0]*V[0,0].conjugate()*inverse(Omega-b1)+V[0,1]*V[0,1].conjugate()*inverse(Omega-b2))/2.0
+    delta[0,1] << (V[0,0]*V[1,0].conjugate()*inverse(Omega-b1)+V[0,1]*V[1,1].conjugate()*inverse(Omega-b2))/2.0
+    delta[1,0] << (V[1,0]*V[0,0].conjugate()*inverse(Omega-b1)+V[1,1]*V[0,1].conjugate()*inverse(Omega-b2))/2.0
+    delta[1,1] << (V[1,0]*V[1,0].conjugate()*inverse(Omega-b1)+V[1,1]*V[1,1].conjugate()*inverse(Omega-b2))/2.0
+    # construct G
+    G[0].zero()
+    for i in range(0,10,2):
+        G[0]['ud'][i:i+2,i:i+2] << inverse(Omega-delta)
+    G[0]['ud'] << inverse(inverse(G[0]['ud']) - Hloc)
+
+    # transform each block using a random transformation matrix
+    for i in range(0,10,2):
+        T = get_random_transformation(2)
+        G[0]['ud'][i:i+2,i:i+2].from_L_G_R(T, G[0]['ud'][i:i+2,i:i+2], T.conjugate().transpose())
+        # 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)
+
+    # analyse the block structure
+    G_new = SK.analyse_block_structure_from_gf(G)
+
+    assert len(SK.deg_shells[0]) == 2, "wrong number of equivalent groups found"
+    assert sorted([len(d) for d in SK.deg_shells[0]]) == [2,3], "wrong number of members in the equivalent groups found"
+    for d in SK.deg_shells[0]:
+        if len(d)==2:
+            assert 'ud_3' in d, "shell ud_3 missing"
+            assert 'ud_4' in d, "shell ud_4 missing"
+        if len(d)==3:
+            assert 'ud_0' in d, "shell ud_0 missing"
+            assert 'ud_1' in d, "shell ud_1 missing"
+            assert 'ud_2' in d, "shell ud_2 missing"
+
+        # the convention is that for every degenerate shell, the transformation
+        # matrix v and the conjugate bool is saved
+        # then,
+        # maybe_conjugate1( v1^dagger G1 v1 ) = maybe_conjugate2( v2^dagger G2 v2 )
+        # therefore, to test, we calculate
+        #   maybe_conjugate( v^dagger G v )
+        # for all degenerate shells and check that they are all equal
+        normalized_gfs = []
+        for key in d:
+            normalized_gf = G_new[0][key].copy()
+            normalized_gf.from_L_G_R(d[key][0].conjugate().transpose(), G_new[0][key], d[key][0])
+            if d[key][1]:
+                conjugate_in_tau(normalized_gf, in_place=True)
+            normalized_gfs.append(normalized_gf)
+        for i in range(len(normalized_gfs)):
+            for j in range(i+1,len(normalized_gfs)):
+                assert_gfs_are_close(normalized_gfs[i], normalized_gfs[j])