[feat] add Pade Sigma analytic continuation and refine tests

python/solid_dmft/dmft_tools/initial_self_energies.py

@@ -498,7 +498,7 @@ def determine_dc_and_initial_sigma(general_params, advanced_params, sum_k,
 
     # Updates the sum_k object with the Matsubara self-energy
     sum_k.put_Sigma([solvers[icrsh].Sigma_freq for icrsh in range(sum_k.n_inequiv_shells)])
-    
+
     # load sigma as first guess in the hartree solver if applicable
     if general_params['solver_type'] == 'hartree':
         # TODO:

python/solid_dmft/postprocessing/maxent_sigma.py

@@ -257,6 +257,8 @@ def main(external_path, iteration=None, continuator_type='inversion_sigmainf', m
     Main function that reads the Matsubara self-energy from h5, analytically continues it,
     writes the results back to the h5 archive and also returns the results.
 
+    Function parallelizes using MPI over impurities and blocks.
+
     Parameters
     ----------
     external_path : string

python/solid_dmft/postprocessing/pade_sigma.py

@@ -0,0 +1,106 @@
+# pyright: reportUnusedExpression=false
+
+import numpy as np
+
+from triqs.utility import mpi
+from h5 import HDFArchive
+from triqs.gf import Gf, MeshReFreq, BlockGf
+
+from solid_dmft.postprocessing.maxent_sigma import _read_h5
+
+def _write_sigma_omega_to_h5(sigma_w, external_path, iteration):
+    """ Writes real-frequency self energy to h5 archive. """
+    h5_internal_path = 'DMFT_results/' + ('last_iter' if iteration is None
+                                          else f'it_{iteration}')
+
+    with HDFArchive(external_path, 'a') as archive:
+        for i, sigma_imp in enumerate(sigma_w):
+            archive[h5_internal_path][f'Sigma_Refreq_{i}'] = sigma_imp
+
+def _run_pade(sigma_iw_list, n_w, w_min, w_max, n_iw, eta):
+    """
+    Run pade in parallel. Call via main function.
+    """
+    mpi.report('Continuing impurities with blocks:')
+
+    imps_blocks = []
+    sigma_iw_flat_list = []
+
+    # create flattened list of self-energies
+    for i, sigma_iw in enumerate(sigma_iw_list):
+        blocks = list(sigma_iw.indices)
+        mpi.report('- Imp {}: {}'.format(i, blocks))
+        for block in blocks:
+            imps_blocks.append((i, block))
+            sigma_iw_flat_list.append(sigma_iw[block])
+
+    sigma_w_flat_list = []
+    wmesh = MeshReFreq(w_min=w_min,w_max=w_max,n_w=n_w)
+    imps_blocks_indices = np.arange(len(imps_blocks))
+    for i in imps_blocks_indices:
+        sigma_w_flat_list.append(Gf(mesh=wmesh, target_shape=sigma_iw_flat_list[i].target_shape))
+
+    # Runs Pade while parallelizing over impurities and blocks
+    for i in mpi.slice_array(imps_blocks_indices):
+        print(f'Rank {mpi.rank} continuing Σ {i}/{len(imps_blocks)}')
+        sigma_w_flat_list[i].set_from_pade(sigma_iw_flat_list[i],n_points=n_iw, freq_offset=eta)
+
+    # sync Pade data
+    for i in imps_blocks_indices:
+        sigma_w_flat_list[i] = mpi.all_reduce(sigma_w_flat_list[i])
+
+    # Create list of BlockGf
+    sigma_w_list = []
+    for i, sigma_iw in enumerate(sigma_iw_list):
+        block_list = []
+        for block in sigma_iw.indices:
+            block_list.append(sigma_w_flat_list.pop(0))
+        sigma_w_list.append(BlockGf(name_list=list(sigma_iw.indices), block_list=block_list, make_copies=True))
+
+    return sigma_w_list
+
+def main(external_path, n_w, w_min, w_max, n_iw, iteration=None, eta=0.0):
+    """
+    Main function that reads the Matsubara self-energy from h5, analytically continues it,
+    writes the results back to the h5 archive and also returns the results.
+
+    Function parallelizes using MPI over impurities and blocks.
+
+    Parameters
+    ----------
+    external_path : string
+        Path to the h5 archive to read from and write to
+    n_w : int
+        number of real frequencies of the final self-energies returned
+    w_min : float
+        Lower end of range where Sigma is being continued.
+    w_max : float
+        Upper end of range where Sigma is being continued.
+    n_iw : int
+        number of Matsubara frequencies to consider for the Pade approximant
+    iteration : int/string
+        Iteration to read from and write to. Default to last_iter
+    eta : float
+        frequency offset within Pade
+
+    Returns
+    -------
+    sigma_w : list of triqs.gf.BlockGf
+        Sigma(omega) per inequivalent shell
+    """
+
+    sigma_iw = None
+    if mpi.is_master_node():
+        sigma_iw, _, _, _ = _read_h5(external_path, iteration)
+    sigma_iw = mpi.bcast(sigma_iw)
+
+    # run pade in parallel
+    sigma_w = _run_pade(sigma_iw, n_w, w_min, w_max, n_iw, eta)
+
+    mpi.report('Writing results to h5 archive now.')
+    if mpi.is_master_node():
+        _write_sigma_omega_to_h5(sigma_w, external_path, iteration)
+    mpi.report('Finished writing Σ(ω) to archive.')
+
+    return sigma_w
+

test/python/CMakeLists.txt

@@ -1,29 +1,29 @@
 # all pytest unittests
-set (all_pytests 
-     test_afm_mapping 
-     test_interaction_hamiltonian 
-     test_manipulate_chemical_potential.py 
-     test_observables.py 
-     test_read_config.py 
-     test_update_dmft_config.py 
+set (all_pytests
+     test_afm_mapping
+     test_interaction_hamiltonian
+     test_manipulate_chemical_potential.py
+     test_observables.py
+     test_read_config.py
+     test_update_dmft_config.py
      test_update_results_h5.py
     )
 
 foreach(test ${all_pytests})
   get_filename_component(test_name ${test} NAME_WE)
   get_filename_component(test_dir ${test} DIRECTORY)
-  
-  add_test(NAME ${test_name} 
-           COMMAND ${TRIQS_PYTHON_EXECUTABLE} -m pytest -vv ${CMAKE_CURRENT_SOURCE_DIR}/${test_dir}/${test_name}.py 
+
+  add_test(NAME ${test_name}
+           COMMAND ${TRIQS_PYTHON_EXECUTABLE} -m pytest -vv ${CMAKE_CURRENT_SOURCE_DIR}/${test_dir}/${test_name}.py
            WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/${test_dir})
-  
+
   set_property(TEST ${test_name} APPEND PROPERTY ENVIRONMENT PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH})
 endforeach()
 
 # ------------------------------#
 
 # all other tests
-set(all_tests 
+set(all_tests
     test_convergence
     test_matheval
     test_plot_correlated_bands
@@ -38,7 +38,7 @@ set(all_tests
 
   # copy reference data for PCB test
   FILE(COPY test_pcb_ref.h5 DESTINATION ${test_dir})
-  
+
   # copy reference data for respack test
   FILE(COPY respack_sfo_data DESTINATION ${test_dir})
 
@@ -52,7 +52,7 @@ endforeach()
 # ------------------------------#
 
 # integration tests
-set (integration_tests 
+set (integration_tests
      svo_hubbardI_basic
      svo_hartree
      lno_hubbardI_mag
@@ -65,17 +65,17 @@ FILE(COPY UIJKL DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
 
 foreach(test ${integration_tests})
   set (test_dir ${CMAKE_CURRENT_BINARY_DIR}/${test})
- 
+
   foreach(file dmft_config.ini inp.h5 ref.h5 test.py)
     FILE(COPY ${test}/${file} DESTINATION ${test_dir})
   endforeach()
-  
-  add_test(NAME ${test} 
+
+  add_test(NAME ${test}
            #COMMAND bash ${test}.sh
            COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_PREFLAGS} ${TRIQS_PYTHON_EXECUTABLE} test.py
            WORKING_DIRECTORY ${test_dir}
           )
-  
+
   set_property(TEST ${test} APPEND PROPERTY ENVIRONMENT PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH})
 endforeach()
 

test/python/test_maxent.py

@@ -3,18 +3,51 @@
 
 from helper import are_iterables_equal
 
-from solid_dmft.postprocessing import maxent_gf_imp, maxent_gf_latt
+from solid_dmft.postprocessing import maxent_gf_imp, maxent_gf_latt, maxent_sigma, pade_sigma
+
+# Note: no unit tests here because lack of MPI support
 
 # Runs maxent on lattice Green function and compares afterwards
-maxent_gf_latt.main('svo_hubbardI_basic/out/inp.h5', sum_spins=True, n_points_maxent=100, n_points_alpha=25, omega_min=-20, omega_max=20)
+mpi.report('#########\nTesting lattice Gf Maxent\n#########')
+maxent_gf_latt.main('svo_hubbardI_basic/out/inp.h5',
+                    sum_spins=True,
+                    n_points_maxent=100,
+                    n_points_alpha=25,
+                    omega_min=-20,
+                    omega_max=20)
 
 # Runs maxent on the impurity Green function
 # No comparison to reference because the spectral function would be too "spiky"
 # because of HubbardI so that numerical differences can dominate the comparison
+mpi.report('#########\nTesting impurity Gf Maxent\n#########')
 maxent_gf_imp.main('svo_hubbardI_basic/out/inp.h5', sum_spins=True,
                    n_points_maxent=50, n_points_alpha=20)
 
 if mpi.is_master_node():
     print('Comparing Alatt_maxent')
-    with HDFArchive('svo_hubbardI_basic/out/inp.h5', 'r')['DMFT_results']['last_iter'] as out, HDFArchive('svo_hubbardI_basic/ref.h5', 'r')['DMFT_results']['last_iter'] as ref:
-        assert are_iterables_equal(out['Alatt_maxent'], ref['Alatt_maxent'])
+    with HDFArchive('svo_hubbardI_basic/out/inp.h5', 'r') as out, HDFArchive('svo_hubbardI_basic/ref.h5', 'r') as ref:
+        assert are_iterables_equal(out['DMFT_results']['last_iter']['Alatt_maxent'], ref['DMFT_results']['last_iter']['Alatt_maxent'])
+
+# Run sigma maxent
+mpi.report('#########\nTesting Sigma Maxent\n#########')
+maxent_sigma.main(external_path='svo_hubbardI_basic/out/inp.h5',
+                  omega_min=-12, omega_max=12,
+                  maxent_error=0.001, iteration=None,
+                  n_points_maxent=50,
+                  n_points_alpha=10,
+                  analyzer='LineFitAnalyzer',
+                  n_points_interp=501,
+                  n_points_final=501,
+                  continuator_type='inversion_dc')[0]
+
+
+# Run sigma pade
+mpi.report('#########\nTesting Sigma Pade\n#########')
+pade_sigma.main(external_path='svo_hubbardI_basic/out/inp.h5',
+                n_w = 4001,
+                w_min=-4.5,
+                w_max=4.5,
+                n_iw=100,
+                eta=0.0
+                )
+