from pytriqs.gf import * 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 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]: 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)): 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) # for testing symm_deg_gf below, we need this # we construct it so that for every group of degenerate blocks of G[0], the # mean of the blocks of G_noisy is equal to G[0] G_noisy = G[0].copy() noise1 = np.random.randn(*delta.target_shape) G_noisy['ud'][:2,:2].data[:,:,:] += noise1 G_noisy['ud'][2:4,2:4].data[:,:,:] -= noise1/2.0 G_noisy['ud'][4:6,4:6].data[:,:,:] -= noise1/2.0 noise2 = np.random.randn(*delta.target_shape) G_noisy['ud'][6:8,6:8].data[:,:,:] += noise2 G_noisy['ud'][8:,8:].data[:,:,:] -= noise2 # for testing backward-compatibility in symm_deg_gf, we need the # un-transformed Green's functions G_pre_transform = G[0].copy() G_noisy_pre_transform = G_noisy.copy() # 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()) 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]: G[0]['ud'][i:i+2,i:i+2] << G[0]['ud'][i:i+2,i:i+2].transpose() 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, 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_pre_transform = SK.block_structure.convert_gf(G_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) 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]: 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)): # here, we use a threshold that is 1 order of magnitude less strict # because of numerics assert_gfs_are_close(normalized_gfs[i], normalized_gfs[j], 1.e-6) # now we check symm_deg_gf # symmetrizing the GF has is has to leave it unchanged 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, 1.e-6) # symmetrizing the noisy GF, which was carefully constructed, # has to give the same result as G_new[0] SK.symm_deg_gf(G_noisy, 0) assert_block_gfs_are_close(G_new[0], G_noisy, 1.e-6) # check backward compatibility of symm_deg_gf # first, construct the old format of the deg shells for ish in range(len(SK.deg_shells)): for gr in range(len(SK.deg_shells[ish])): SK.deg_shells[ish][gr] = SK.deg_shells[ish][gr].keys() # symmetrizing the GF as is has to leave it unchanged G_new_symm << G_pre_transform SK.symm_deg_gf(G_new_symm, 0) assert_block_gfs_are_close(G_new_symm, G_pre_transform, 1.e-6) # symmetrizing the noisy GF pre transform, which was carefully constructed, # has to give the same result as G_pre_transform SK.symm_deg_gf(G_noisy_pre_transform, 0) assert_block_gfs_are_close(G_noisy_pre_transform, G_pre_transform, 1.e-6)