From a87276b6c4d5994ce6c347dcd6c04699859daa17 Mon Sep 17 00:00:00 2001 From: Alexander Hampel Date: Sun, 26 May 2024 20:35:41 -0400 Subject: [PATCH] [feat] add option to use DLR mesh in Sumk (#254) allow using the MeshDLRImFreq to be used as general Sumk mesh during the DMFT loop for all functions. --- doc/ChangeLog.md | 2 + python/triqs_dft_tools/sumk_dft.py | 26 ++++++++----- test/python/calc_mu.py | 8 +++- test/python/srvo3_Gloc.py | 59 ++++++++++++++++++++++++++---- 4 files changed, 76 insertions(+), 19 deletions(-) diff --git a/doc/ChangeLog.md b/doc/ChangeLog.md index a42fbaf8..ac2f2390 100644 --- a/doc/ChangeLog.md +++ b/doc/ChangeLog.md @@ -8,6 +8,7 @@ DFTTools Version 3.3.0 is a release that * is compatible with TRIQS 3.3.x * includes the latest app4triqs changes * introduce `dc_imp_dyn` attribute in sumk object to store dynamic part of DC potential +* allows using MeshDLRImFreq as Sumk mesh * improved standard behavior of block struct (#248) (see below for details) We thank all contributors: Sophie Beck, Thomas Hahn, Alexander Hampel, Henri Menke, Dylan Simon, Nils Wentzell @@ -21,6 +22,7 @@ Find below an itemized list of changes in this release. ### feat * allow dict/np.ndarrays input in `symm_deg_gf` * introduce `dc_imp_dyn` attribute in sumk object to store dynamic part of DC potential +* allows using MeshDLRImFreq as Sumk mesh * previously the default `gf_struct_solver` in a initialized blockstructure had keys `up` / `down`, inconsistent with the default behavior after running `analyse_block_structure`: `up_0` / `down_0`. Now the default solver structure always has the `_0` in the key. * old behavior resulted in error when analyse_block_structure was called diff --git a/python/triqs_dft_tools/sumk_dft.py b/python/triqs_dft_tools/sumk_dft.py index 9aeceae4..c5550b16 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..09354211 100644 --- a/test/python/srvo3_Gloc.py +++ b/test/python/srvo3_Gloc.py @@ -19,19 +19,26 @@ # ################################################################################ -from h5 import * -from triqs.gf import * -from triqs_dft_tools.sumk_dft import * -from triqs_dft_tools.converters.wien2k import * +from h5 import HDFArchive +from triqs.utility import mpi +from triqs.gf import MeshImFreq, MeshDLRImFreq, Gf, BlockGf, make_gf_dlr, make_gf_imfreq +from triqs_dft_tools.sumk_dft import SumkDFT from triqs.operators.util import set_operator_structure -from triqs.utility.comparison_tests import * +from triqs.utility.comparison_tests import assert_block_gfs_are_close 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")