save meshes in npz

src/dcore/dcore.py

@@ -57,23 +57,6 @@ def dcore(filename, np=1):
 
     solver.do_steps(max_step=params["control"]["max_step"])
 
-    # Write information for analytic continuation of the self-energy
-    with open(params["model"]["seedname"] + "_anacont.toml", "w") as f:
-        toml.dump(
-            {
-                "beta":       params["system"]["beta"],
-                "Nomega":     params["tool"]["Nomega"],
-                "omega_min":  params["tool"]["omega_min"],
-                "omega_max":  params["tool"]["omega_max"],
-                "pade" : {
-                    "n_min"    : params["tool"]["n_pade_min"],
-                    "n_max"    : params["tool"]["n_pade_max"],
-                    "omega_max": params["tool"]["omega_pade"],
-                    "eta"      : params["tool"]["eta"],
-                }
-            },
-        f)
-
     print("\n########################  Done  ########################\n")
 
 

src/dcore/dcore_anacont.py

@@ -27,17 +27,17 @@
 from dcore.anacont import pade, spm
 
 def dcore_anacont(inifile):
-    parser = create_parser(["system", "model", "impurity_solver", "post", "post.anacont"])
+    # tool is removed but read for error message
+    parser = create_parser(["system", "model", "impurity_solver", "tool", "post", "post.anacont"])
     parser.read(inifile)
     params = parser.as_dict()
     parse_parameters(params)
 
     beta = params["system"]["beta"]
 
-    omega_min = params["post"]["omega_min"]
-    omega_max = params["post"]["omega_max"]
-
-    Nomega = params["post"]["Nomega"]
+    omega_min = params["post.anacont"]["omega_min"]
+    omega_max = params["post.anacont"]["omega_max"]
+    Nomega = params["post.anacont"]["Nomega"]
     mesh_w = MeshReFreq(omega_min, omega_max, Nomega)
 
     seedname = params["model"]["seedname"]
@@ -60,11 +60,15 @@ def dcore_anacont(inifile):
         data_w = spm.anacont(sigma_iw_npz, beta, mesh_w, params_ac)
     else:
         assert False, "Unknown solver: " + solver
-
+
+    data_w["omega"] = mesh_w.values()
     if params["post.anacont"]["show_result"] or params["post.anacont"]["save_result"]:
         import matplotlib.pyplot as plt
 
-        ndata = len(data_w)
+        ndata = 0
+        for key in data_w.keys():
+            if key.startswith("data"):
+                ndata += 1
         for idata in range(ndata):
             sigma_w = data_w[f"data{idata}"]
             for iorb, jorb in itertools.product(range(sigma_w.shape[1]), range(sigma_w.shape[2])):

src/dcore/dcore_check.py

@@ -39,8 +39,9 @@ def __init__(self, ini_file, max_n_iter):
         print("\n  @ Reading {0} ...".format(ini_file))
         #
         # Construct a parser with default values
+        # tool is removed but read for error message
         #
-        pars = create_parser(['model', 'system', 'tool'])
+        pars = create_parser(['model', 'system', 'tool', 'post.check'])
         #
         # Parse keywords and store
         #
@@ -51,7 +52,6 @@ def __init__(self, ini_file, max_n_iter):
         # Delete unnecessary parameters
         delete_parameters(self.p, block='model', retain=['seedname'])
         delete_parameters(self.p, block='system', retain=['beta', 'n_iw', 'mu', 'fix_mu'])
-        delete_parameters(self.p, block='tool', retain=['omega_check'])
 
         # Summary of input parameters
         print_parameters(self.p)
@@ -74,7 +74,7 @@ def __init__(self, ini_file, max_n_iter):
         print("  Total number of Iteration: {0}".format(self.n_iter))
 
         # If omega_check is not specified, a fixed number of Matsubara points are taken
-        self.omega_check = self.p['tool']['omega_check']
+        self.omega_check = self.p['post.check']['omega_check']
         if self.omega_check == 0:
             nmax = min(30, self.n_iw)
             self.omega_check = (2*nmax+1) * math.pi / self.beta

src/dcore/dcore_spectrum.py

@@ -37,6 +37,10 @@
 from dcore.lattice_models import create_lattice_model
 from .sumkdft_workers.launcher import run_sumkdft
 
+def _read_ws(npz_file):
+    npz = np.load(npz_file)
+    ws = npz["omega"]
+    return ws
 
 def _read_sigma_w(npz_file, nsh, mesh, block_names):
     npz = np.load(npz_file)
@@ -160,10 +164,18 @@ def __init__(self, seedname, params, n_k, xk, nkdiv_mesh, kvec_mesh, dir_post):
         """
 
         self._params = copy.deepcopy(params)
+        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
+
         # Construct a SumKDFT object
-        self._omega_min = float(params["post"]["omega_min"])
-        self._omega_max = float(params["post"]["omega_max"])
-        self._Nomega = int(params["post"]["Nomega"])
+        self._omega_min = ws[0]
+        self._omega_max = ws[-1]
+        self._Nomega = len(ws)
+        # self._omega_min = float(params["post"]["omega_min"])
+        # self._omega_max = float(params["post"]["omega_max"])
+        # self._Nomega = int(params["post"]["Nomega"])
         self._broadening = float(params["post.spectrum"]["broadening"])
         self._seedname = seedname
         self._n_k = n_k
@@ -172,7 +184,6 @@ def __init__(self, seedname, params, n_k, xk, nkdiv_mesh, kvec_mesh, dir_post):
         self._nkdiv_mesh = nkdiv_mesh
         self._dir_post = dir_post
 
-        self._sigma_w_file = os.path.join(dir_post, "sigma_w.npz")
 
         self._solver = DMFTPostSolver(
             seedname, self._params, output_file=seedname + ".out.h5"
@@ -269,13 +280,12 @@ def post(self):
         #
         # Real-frequency self-energy
         #
-        mesh = [self._omega_min, self._omega_max, self._Nomega]
         if not os.path.exists(self._sigma_w_file):
             sys.exit("File not found: " + self._sigma_w_file)
         sigma_w_sh = _read_sigma_w(
             self._sigma_w_file,
             self._solver.n_inequiv_shells,
-            MeshReFreq(*mesh),
+            self.mesh_w,
             self._solver.spin_block_names,
         )
 
@@ -286,6 +296,7 @@ def post(self):
         print("\n Writing real-freqnecy self-energy into ", filename)
         save_Sigma_w_sh_txt(filename, sigma_w_sh, self._solver.spin_block_names)
 
+        mesh = [self._omega_min, self._omega_max, self._Nomega]
         #
         #  (Partial) DOS
         #
@@ -474,6 +485,7 @@ def dcore_spectrum(filename, np=1):
     print("  Input File Name : ", filename)
     #
     # Construct a parser with default values
+    # tool is removed but read for error message
     #
     pars = create_parser(["model", "system", "impurity_solver", "tool", "post", "post.spectrum", "mpi"])
     #
@@ -512,8 +524,6 @@ def dcore_spectrum(filename, np=1):
     # make directory
     dir_post = p["post"]["dir_post"]
 
-    # for backward compatibility
-    p["tool"]["post_dir"] = dir_post
 
     if dir_post:
         os.makedirs(dir_post, exist_ok=True)

src/dcore/dmft_core.py

@@ -943,6 +943,9 @@ def _save_sigma_iw(self, dm_sh: List[numpy.ndarray]) -> None:
                 data[f"data{idata}"] = g.data
                 data[f"hartree_fock{idata}"] = hf[bname]
                 idata += 1
+                iwn = g.mesh.values()
+        data["beta"] = self._beta
+        data["iwn"] = iwn
         numpy.savez(self._seedname + "_sigma_iw.npz", **data)
 
     def chemical_potential(self, iteration_number):