From 8d6d8b53c5b73c1033a6f7b5e742b650982a7981 Mon Sep 17 00:00:00 2001 From: "Gernot J. Kraberger" Date: Mon, 19 Mar 2018 11:09:31 +0100 Subject: [PATCH] SumkDFT: analyze_block_structure_from_gf for Gf Re/Im Time/Freq --- python/sumk_dft.py | 114 ++++++++++++---------- test/CMakeLists.txt | 2 +- test/analyze_block_structure_from_gf.py | 10 +- test/analyze_block_structure_from_gf2.py | 115 +++++++++++++++++++++++ 4 files changed, 187 insertions(+), 54 deletions(-) create mode 100644 test/analyze_block_structure_from_gf2.py diff --git a/python/sumk_dft.py b/python/sumk_dft.py index 33fb8b23..7f527441 100644 --- a/python/sumk_dft.py +++ b/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: - assert all(isinstance(g_sh._first(), GfImTime) for g_sh in G), "G must be a BlockGf of either GfImFreq or GfImTime" + # keep a GfImTime from the supplied 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: - assert all(isinstance(g_sh._first(), GfImTime) for g_sh in G), "G must be a BlockGf of either GfImFreq or GfImTime" + # keep a GfImTime from the supplied 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, 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) + gf_to_symm[key].from_L_G_R(v, ss.transpose(), v.conjugate().transpose()) 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 diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 02edb21a..b1dc6e8c 100644 --- a/test/CMakeLists.txt +++ b/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) diff --git a/test/analyze_block_structure_from_gf.py b/test/analyze_block_structure_from_gf.py index 94ec6e18..c777e9c3 100644 --- a/test/analyze_block_structure_from_gf.py +++ b/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) - conjugate_in_tau(G_noisy['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() + 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)): diff --git a/test/analyze_block_structure_from_gf2.py b/test/analyze_block_structure_from_gf2.py new file mode 100644 index 00000000..d371a9c5 --- /dev/null +++ b/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"