[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

@@ -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

python/triqs_dft_tools/sumk_dft.py

@@ -58,7 +58,7 @@ def __init__(self, hdf_file, h_field=0.0, mesh=None, beta=40, n_iw=1025, use_dft
                   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 @@ def __init__(self, hdf_file, h_field=0.0, mesh=None, beta=40, n_iw=1025, use_dft
                 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,16 +566,20 @@ def lattice_gf(self, ik, mu=None, broadening=None, mesh=None, with_Sigma=True, w
                 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:
                 mesh = self.mesh
                 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 @@ def lattice_gf(self, ik, mu=None, broadening=None, mesh=None, with_Sigma=True, w
             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 @@ def lattice_gf(self, ik, mu=None, broadening=None, mesh=None, with_Sigma=True, w
         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 @@ def put_Sigma(self, Sigma_imp, transform_to_sumk_blocks=True):
         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)]

test/python/calc_mu.py

@@ -31,7 +31,8 @@
 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 @@ def test_dichotomy(self):
         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)

test/python/srvo3_Gloc.py

@@ -19,35 +19,78 @@
 #
 ################################################################################
 
-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
+
+# classic full Matsubara mesh
+mpi.report(f"{'#'*12}\nregular Matsubara mesh test\n")
 
-# Init the SumK class
-SK=SumkDFT(hdf_file='SrVO3.ref.h5',use_dft_blocks=True)
+# 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']
 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")