make dcore_anacont and spectrum work with TRIQS

src/dcore/anacont/spm.py

@@ -334,9 +334,11 @@ def anacont(sigma_iw_npz, beta, mesh_w, params_spm, params_spm_solver):
         mesh_iw = MeshImFreq(beta, "Fermion", n_matsubara_retain)
         gf_iw = GfImFreq(data=data, beta=beta, mesh=mesh_iw)
         n_orbitals = gf_iw.data.shape[1]
-        matsubara_frequencies = np.imag(gf_iw.mesh.values()[n_matsubara_retain:])
+        iws = list(gf_iw.mesh.values())
+        ws = list(mesh_w.values())
+        matsubara_frequencies = np.imag(iws[n_matsubara_retain:])
         sigma_w_data = np.zeros(
-            (mesh_w.size, n_orbitals, n_orbitals), dtype=np.complex128
+            (len(ws), n_orbitals, n_orbitals), dtype=np.complex128
         )
         for i_orb in range(n_orbitals):
             print(
@@ -351,7 +353,7 @@ def anacont(sigma_iw_npz, beta, mesh_w, params_spm, params_spm_solver):
             )
 
             density, gf_tau_fit, energies_extract, sum_rule_const, chi2 = get_single_continuation(
-                tau_grid, gf_tau, n_sv, beta, mesh_w.values()[0], mesh_w.values()[-1], mesh_w.size, sum_rule_const=const_imag_tail, lambd=L1_coeff, verbose=verbose_opt, solver=solver, solver_opts=params_spm_solver
+                tau_grid, gf_tau, n_sv, beta, ws[0], ws[-1], len(ws), sum_rule_const=const_imag_tail, lambd=L1_coeff, verbose=verbose_opt, solver=solver, solver_opts=params_spm_solver
             )
             energies, gf_real, gf_imag = get_kramers_kronig_realpart(
                 energies_extract, dos_to_gf_imag(density)
@@ -373,4 +375,4 @@ def anacont(sigma_iw_npz, beta, mesh_w, params_spm, params_spm_solver):
                 plt.close()
         sigma_w = GfReFreq(mesh=mesh_w, data=sigma_w_data)
         data_w[key] = sigma_w.data
-    return data_w
+    return data_w

src/dcore/dcore_anacont.py

@@ -61,9 +61,17 @@ def dcore_anacont(inifile):
             params_ac["post.anacont.spm.solver"] = {}
         data_w = spm.anacont(sigma_iw_npz, beta, mesh_w, params_ac["post.anacont.spm"], params_ac["post.anacont.spm.solver"])
     else:
-        assert False, "Unknown solver: " + solver
+        sys.exit("ERROR: Unknown anacont solver " + solver)
+    omega = numpy.array(list(mesh_w.values()))
+    data_w["omega"] = omega
+
+    dir_post = params["post"]["dir_post"]
+    if not os.path.exists(dir_post):
+        os.makedirs(dir_post)
+    file_sigma_w = os.path.join(dir_post, "sigma_w.npz")
+    print("Writing to", file_sigma_w + "...")
+    numpy.savez(file_sigma_w, **data_w)
 
-    data_w["omega"] = mesh_w.values()
     if params["post.anacont"]["show_result"] or params["post.anacont"]["save_result"]:
         import matplotlib.pyplot as plt
 
@@ -75,8 +83,8 @@ def dcore_anacont(inifile):
             sigma_w = data_w[f"data{idata}"]
             for iorb, jorb in itertools.product(range(sigma_w.shape[1]), range(sigma_w.shape[2])):
                 plt.axhline(y=0, xmin=omega_min, xmax=omega_max, color="lightgrey")
-                plt.plot(mesh_w.points, sigma_w[:,iorb,jorb].real, label="Real")
-                plt.plot(mesh_w.points, sigma_w[:,iorb,jorb].imag, label="Imag")
+                plt.plot(omega, sigma_w[:,iorb,jorb].real, label="Real")
+                plt.plot(omega, sigma_w[:,iorb,jorb].imag, label="Imag")
                 plt.xlim(omega_min, omega_max)
                 plt.xlabel(r"$\omega$")
                 plt.legend()
@@ -90,12 +98,6 @@ def dcore_anacont(inifile):
                     plt.show()
                 plt.close()
 
-    dir_post = params["post"]["dir_post"]
-    if not os.path.exists(dir_post):
-        os.makedirs(dir_post)
-    file_sigma_w = os.path.join(dir_post, "sigma_w.npz")
-    print("Writing to", file_sigma_w + "...")
-    numpy.savez(file_sigma_w, **data_w)
 
 def run():
 
@@ -120,4 +122,4 @@ def run():
 
     args = parser.parse_args()
 
-    dcore_anacont(args.path_input_files)
+    dcore_anacont(args.path_input_files)

src/dcore/dcore_spectrum.py

@@ -167,7 +167,12 @@ def __init__(self, seedname, params, n_k, xk, nkdiv_mesh, kvec_mesh, dir_post):
         self._sigma_w_file = os.path.join(dir_post, "sigma_w.npz")
         ws = _read_ws(self._sigma_w_file)
         self.mesh_w = MeshReFreq(ws[0], ws[-1], len(ws))
-        self.mesh_w._points = ws
+        if TRIQS_COMPAT:
+            self.mesh_w._points = ws
+        else:
+            ws_linear = np.linspace(ws[0], ws[-1], num=len(ws))
+            if not np.allclose(ws, ws_linear):
+                sys.exit("Error: omega points must be linearly spaced when using TRIQS (DCORE_TRIQS_COMPAT=0)")
 
         # Construct a SumKDFT object
         self._omega_min = ws[0]