diff --git a/python/triqs_dft_tools/sumk_dft.py b/python/triqs_dft_tools/sumk_dft.py index 676ad395..1f046070 100644 --- a/python/triqs_dft_tools/sumk_dft.py +++ b/python/triqs_dft_tools/sumk_dft.py @@ -58,7 +58,7 @@ class SumkDFT(object): The value of magnetic field to add to the DFT Hamiltonian. The contribution -h_field*sigma is added to diagonal elements of the Hamiltonian. It cannot be used with the spin-orbit coupling on; namely h_field is set to 0 if self.SO=True. - mesh: MeshImFreq or MeshReFreq, optional. Frequency mesh of Sigma. + mesh: MeshImFreq, MeshDLRImFreq, or MeshReFreq, optional. Frequency mesh of Sigma. beta : real, optional Inverse temperature. Used to construct imaginary frequency if mesh is not given. n_iw : integer, optional @@ -115,6 +115,9 @@ class SumkDFT(object): self.mesh_values = np.linspace(self.mesh(self.mesh.first_index()), self.mesh(self.mesh.last_index()), len(self.mesh)) + elif isinstance(mesh, MeshDLRImFreq): + self.mesh = mesh + self.mesh_values = np.array([iwn.value for iwn in mesh.values()]) elif isinstance(mesh, MeshReFreq): self.mesh = mesh self.mesh_values = np.linspace(self.mesh.w_min, self.mesh.w_max, len(self.mesh)) @@ -563,6 +566,8 @@ class SumkDFT(object): mesh = Sigma_imp[0].mesh if isinstance(mesh, MeshImFreq): mesh_values = np.linspace(mesh(mesh.first_index()), mesh(mesh.last_index()), len(mesh)) + elif isinstance(mesh, MeshDLRImFreq): + mesh_values = np.array([iwn.value for iwn in mesh.values()]) else: mesh_values = np.linspace(mesh.w_min, mesh.w_max, len(mesh)) else: @@ -570,9 +575,11 @@ class SumkDFT(object): mesh_values = self.mesh_values elif not mesh is None: - assert isinstance(mesh, MeshReFreq) or isinstance(mesh, MeshImFreq), "mesh must be a triqs MeshReFreq or MeshImFreq" + assert isinstance(mesh, (MeshReFreq, MeshDLRImFreq, MeshImFreq)), "mesh must be a triqs MeshReFreq or MeshImFreq" if isinstance(mesh, MeshImFreq): mesh_values = np.linspace(mesh(mesh.first_index()), mesh(mesh.last_index()), len(mesh)) + elif isinstance(mesh, MeshDLRImFreq): + mesh_values = np.array([iwn.value for iwn in mesh.values()]) else: mesh_values = np.linspace(mesh.w_min, mesh.w_max, len(mesh)) else: @@ -586,12 +593,8 @@ class SumkDFT(object): gf_struct = [(spn[isp], block_structure[isp]) for isp in range(self.n_spin_blocks[self.SO])] block_ind_list = [block for block, inner in gf_struct] - if isinstance(mesh, MeshImFreq): - glist = lambda: [Gf(mesh=mesh, target_shape=[len(inner),len(inner)]) - for block, inner in gf_struct] - else: - glist = lambda: [Gf(mesh=mesh, target_shape=[len(inner),len(inner)]) - for block, inner in gf_struct] + glist = lambda: [Gf(mesh=mesh, target_shape=[len(inner),len(inner)]) + for block, inner in gf_struct] G_latt = BlockGf(name_list=block_ind_list, block_list=glist(), make_copies=False) G_latt.zero() @@ -603,7 +606,7 @@ class SumkDFT(object): for ibl, (block, gf) in enumerate(G_latt): ind = ntoi[spn[ibl]] n_orb = self.n_orbitals[ik, ind] - if isinstance(mesh, MeshImFreq): + if isinstance(mesh, (MeshImFreq, MeshDLRImFreq)): gf.data[:, :, :] = (idmat[ibl] * (mesh_values[:, None, None] + mu + self.h_field*(1-2*ibl)) - self.hopping[ik, ind, 0:n_orb, 0:n_orb]) else: @@ -647,7 +650,10 @@ class SumkDFT(object): assert len(Sigma_imp) == self.n_corr_shells,\ "put_Sigma: give exactly one Sigma for each corr. shell!" - if isinstance(self.mesh, MeshImFreq) and all(isinstance(gf.mesh, MeshImFreq) and isinstance(gf, Gf) and gf.mesh == self.mesh for bname, gf in Sigma_imp[0]): + if (isinstance(self.mesh, (MeshImFreq, MeshDLRImFreq)) and + all(isinstance(gf.mesh, (MeshImFreq, MeshDLRImFreq)) and + isinstance(gf, Gf) and + gf.mesh == self.mesh for bname, gf in Sigma_imp[0])): # Imaginary frequency Sigma: self.Sigma_imp = [self.block_structure.create_gf(ish=icrsh, mesh=Sigma_imp[icrsh].mesh, space='sumk') for icrsh in range(self.n_corr_shells)] diff --git a/test/python/calc_mu.py b/test/python/calc_mu.py index ae98fcca..b2c90961 100644 --- a/test/python/calc_mu.py +++ b/test/python/calc_mu.py @@ -31,7 +31,8 @@ from triqs_dft_tools.sumk_dft import * class test_solver(unittest.TestCase): def setUp(self): - self.iw_mesh = MeshImFreq(beta=40, S='Fermion', n_iw=300) + self.iw_mesh = MeshImFreq(beta=40, statistic='Fermion', n_iw=300) + self.dlr_mesh = MeshDLRImFreq(beta=40, statistic='Fermion', w_max=10, eps=1e-10) self.w_mesh = MeshReFreq(n_w=1001, window=(-3,3)) # magic reference value for the Wien2k SVO t2g example self.ref_mu = 0.281 @@ -42,6 +43,11 @@ class test_solver(unittest.TestCase): mu = sumk.calc_mu(method='dichotomy', precision=0.001, delta=0.1) self.assertTrue(abs(self.ref_mu - mu) < 0.01) + def test_dichotomy_dlr(self): + sumk = SumkDFT('SrVO3.ref.h5', mesh=self.dlr_mesh) + mu = sumk.calc_mu(method='dichotomy', precision=0.001, delta=0.1) + self.assertTrue(abs(self.ref_mu - mu) < 0.01) + def test_dichotomy_real(self): sumk = SumkDFT('SrVO3.ref.h5', mesh=self.w_mesh) mu = sumk.calc_mu(method='dichotomy', precision=0.001, delta=0.1, broadening = 0.01, beta=1000) diff --git a/test/python/srvo3_Gloc.py b/test/python/srvo3_Gloc.py index ebeaf479..083b6f9b 100644 --- a/test/python/srvo3_Gloc.py +++ b/test/python/srvo3_Gloc.py @@ -20,18 +20,25 @@ ################################################################################ from h5 import * -from triqs.gf import * +from triqs.gf import MeshImFreq, MeshDLRImFreq from triqs_dft_tools.sumk_dft import * from triqs_dft_tools.converters.wien2k import * from triqs.operators.util import set_operator_structure from triqs.utility.comparison_tests import * from triqs.utility.h5diff import h5diff +import time + # Basic input parameters beta = 40 +n_iw = 1025 -# Init the SumK class -SK=SumkDFT(hdf_file='SrVO3.ref.h5',use_dft_blocks=True) +# classic full Matsubara mesh +mpi.report(f"{'#'*12}\nregular Matsubara mesh test\n") + +# Init the SumK class (reference data with n_iw=1025) +iw_mesh = MeshImFreq(n_iw=n_iw,beta=beta, statistic='Fermion') +SK=SumkDFT(hdf_file='SrVO3.ref.h5',mesh=iw_mesh,use_dft_blocks=True) num_orbitals = SK.corr_shells[0]['dim'] l = SK.corr_shells[0]['l'] @@ -39,15 +46,51 @@ spin_names = ['down','up'] orb_hybridized = False gf_struct = set_operator_structure(spin_names,num_orbitals,orb_hybridized) -glist = [ GfImFreq(target_shape=(bl_size,bl_size),beta=beta) for bl, bl_size in gf_struct] +glist = [ Gf(target_shape=(bl_size,bl_size),mesh=iw_mesh) for bl, bl_size in gf_struct] Sigma_iw = BlockGf(name_list = [bl for bl, bl_size in gf_struct], block_list = glist, make_copies = False) SK.set_Sigma([Sigma_iw]) + +if mpi.is_master_node(): + start_time = time.time() + Gloc = SK.extract_G_loc() +if mpi.is_master_node(): + mpi.report(f'extract_G_loc time: {(time.time()-start_time)*1000:.1f} msec') + if mpi.is_master_node(): with HDFArchive('srvo3_Gloc.out.h5','w') as ar: ar['Gloc'] = Gloc[0] if mpi.is_master_node(): h5diff("srvo3_Gloc.out.h5","srvo3_Gloc.ref.h5") + +mpi.report(f"{'#'*12}\n") + + +# DLR Matsubara mesh +mpi.report(f"{'#'*12}\nDLR Matsubara mesh test\n") + +dlr_mesh = MeshDLRImFreq(beta=beta, statistic='Fermion', w_max=10, eps=1e-10) +SK=SumkDFT(hdf_file='SrVO3.ref.h5',mesh=dlr_mesh,use_dft_blocks=True) + +glist_dlr = [ Gf(target_shape=(bl_size,bl_size),mesh=dlr_mesh) for bl, bl_size in gf_struct] +Sigma_dlr = BlockGf(name_list = [bl for bl, bl_size in gf_struct], block_list = glist_dlr, make_copies = False) +SK.set_Sigma([Sigma_dlr]) + +if mpi.is_master_node(): + start_time = time.time() + +Gloc_dlr_iw = SK.extract_G_loc() + +if mpi.is_master_node(): + mpi.report(f'extract_G_loc time: {(time.time()-start_time)*1000:.1f} msec') + + with HDFArchive('srvo3_Gloc.out.h5','a') as ar: + ar['Gloc_dlr'] = make_gf_imfreq(make_gf_dlr(Gloc_dlr_iw[0]),n_iw=n_iw) + # get full Giw and compare + Gloc_iw_full = make_gf_imfreq(make_gf_dlr(Gloc_dlr_iw[0]),n_iw=n_iw) + assert_block_gfs_are_close(Gloc[0], Gloc_iw_full) + +mpi.report(f"{'#'*12}\n")