diff --git a/doc/reference/input.rst b/doc/reference/input.rst
index 6c8c2ced..d64c5a92 100644
--- a/doc/reference/input.rst
+++ b/doc/reference/input.rst
@@ -16,13 +16,14 @@ The following table shows which blocks are used by each program.
     :widths: 5, 5, 5, 5, 5
 
     [model]          ,  Yes, Yes, Yes, Yes
+    [pre]            ,  Yes,    ,    ,
     [system]         ,     , Yes, Yes, Yes
     [impurity_solver],     , Yes,    , Yes
     [control]        ,     , Yes
     [tool]           ,     ,    , Yes, Yes
     [mpi]            ,     , Yes,    , Yes
 
-For example, we can see that ``dcore_pre`` reads only [model] block. Therefore, ``dcore_pre`` needs to be re-executed only when the parameters in [model] block are changed.
+For example, we can see that ``dcore_pre`` reads only [model] and [pre] blocks. Therefore, ``dcore_pre`` needs to be re-executed only when the parameters in [model] and [pre] blocks are changed.
 
 The parameters included in each block are explained below.
 
@@ -54,6 +55,13 @@ See separate pages from the link below for detailed descriptions of some paramet
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 
+[pre] block
+--------------
+
+This block includes parameters that controls the behavior of ``dcore_pre``.
+
+.. include:: pre_desc.txt
+
 [system] block
 --------------
 
diff --git a/src/dcore/dcore_bse.py b/src/dcore/dcore_bse.py
index a9467aa9..531a8959 100644
--- a/src/dcore/dcore_bse.py
+++ b/src/dcore/dcore_bse.py
@@ -27,7 +27,7 @@
 
 from dcore._dispatcher import HDFArchive, dyson
 from dcore.dmft_core import DMFTCoreSolver
-from dcore.program_options import create_parser, parse_parameters
+from dcore.program_options import create_parser, parse_parameters, print_parameters, delete_parameters
 from dcore.tools import *
 from dcore import impurity_solvers
 from .sumkdft_workers.launcher import run_sumkdft
@@ -42,17 +42,14 @@ def to_str(x):
     return x
 
 
-def compare_str_list(list1, list2):
-    if len(list1) != len(list2):
-        return False
+def _compare_str_list(list1, list2):
+    assert len(list1) == len(list2), f"len({list1}) != len({list2})"
     for x, y in zip(list1, list2):
-        if to_str(x) != to_str(y):
-            return False
-    return True
+        assert to_str(x) == to_str(y), f"{to_str(x)} != {to_str(y)}"
 
 
 def calc_g2_in_impurity_model(solver_name, solver_params, mpirun_command, basis_rot, Umat, gf_struct, beta, n_iw,
-                              Sigma_iw, Gloc_iw, num_wb, num_wf, ish, freqs=None):
+                              Sigma_iw, Gloc_iw, num_wb, num_wf, only_chiloc, ish, freqs=None):
     """
 
     Calculate G2 in an impurity model
@@ -86,24 +83,29 @@ def calc_g2_in_impurity_model(solver_name, solver_params, mpirun_command, basis_
     # Solve the model
     rot = impurity_solvers.compute_basis_rot(basis_rot, sol)
     if flag_box:
-        xloc, chiloc = sol.calc_Xloc_ph(rot, mpirun_command, num_wf, num_wb, s_params)
+        xloc, chiloc = sol.calc_Xloc_ph(rot, mpirun_command, num_wf, num_wb, s_params, only_chiloc)
     else:
         xloc, chiloc = sol.calc_Xloc_ph_sparse(rot, mpirun_command, freqs, num_wb, s_params)
 
     # Check results for x_loc
     print("\n Checking x_loc...")
-    assert isinstance(xloc, dict)
-    for key, data in list(xloc.items()):
-        # print("  ", key)
-        if flag_box:
-            assert data.shape == (num_wb, 2*num_wf, 2*num_wf)
-        else:
-            assert data.shape == (freqs.shape[0],)
-    print(" OK")
+    if xloc is None:
+        print("   not computed")
+    else:
+        assert isinstance(xloc, dict)
+        for key, data in list(xloc.items()):
+            # print("  ", key)
+            if flag_box:
+                assert data.shape == (num_wb, 2*num_wf, 2*num_wf)
+            else:
+                assert data.shape == (freqs.shape[0],)
+        print(" OK")
 
     # Check results for chi_loc
-    if chiloc is not None:
-        print("\n Checking chi_loc...")
+    print("\n Checking chi_loc...")
+    if chiloc is None:
+        print("   not computed")
+    else:
         assert isinstance(chiloc, dict)
         for key, data in list(chiloc.items()):
             # print("  ", key)
@@ -274,7 +276,7 @@ def __init__(self, h5_file, bse_info, n_corr_shells, n_flavors, use_spin_orbit,
         assert n_block == len(spin_names)
 
         # NOTE: change the order of spins in HDF5 to meet SumkDFTChi
-        self.block2 = IndexPair2(list(range(n_corr_shells)), sorted(spin_names), only_diagonal1=only_diagonal)
+        self.block2 = IndexPair2(list(range(n_corr_shells)), spin_names, only_diagonal1=only_diagonal)
         self.inner2 = IndexPair(list(range(n_inner)), convert_to_int=True)
         print(" block2 namelist =", self.block2.namelist)
         print(" inner2 namelist =", self.inner2.namelist)
@@ -284,8 +286,8 @@ def __init__(self, h5_file, bse_info, n_corr_shells, n_flavors, use_spin_orbit,
         self.h5bse = h5BSE(h5_file, bse_grp)
         if bse_info == 'check':
             # check equivalence of info
-            #assert compare_str_list(h5bse.get(key=('block_name',)), self.block2.namelist)
-            #assert compare_str_list(h5bse.get(key=('inner_name',)), self.inner2.namelist)
+            _compare_str_list(self.h5bse.get(key=('block_name',)), self.block2.namelist)
+            _compare_str_list(self.h5bse.get(key=('inner_name',)), self.inner2.namelist)
             assert self.h5bse.get(key=('beta',)) == beta
         elif bse_info == 'save':
             # save info
@@ -337,7 +339,7 @@ def save_xloc(self, xloc_ijkl, icrsh):
         self._save_common(xloc_ijkl, icrsh, 'X_loc')
 
     def save_chiloc(self, chiloc_ijkl, icrsh):
-        self._save_common(chiloc_ijkl, icrsh, 'chi_loc')
+        self._save_common(chiloc_ijkl, icrsh, 'chi_loc_in')
 
     def save_gamma0(self, u_mat, icrsh):
 
@@ -406,7 +408,7 @@ def _calc_bse_x0q(self):
         Calc X_0(q)
         """
         print("\n--- dcore_bse - X_0(q)")
-        if self._params['bse']['skip_X0q_if_exists'] and os.path.exists(self._params['bse']['h5_output_file']):
+        if self._params['bse']['skip_X0q']:
             print(" skip")
             return
 
@@ -452,10 +454,16 @@ def calc_num_flavors(_ish):
         # init for saving data into HDF5
         #
         print("\n--- dcore_bse - invoking h5BSE...")
+
+        if os.path.isfile(self._params['bse']['h5_output_file']):
+            bse_info = 'check'
+        else:
+            bse_info = 'save'
+
         bse = SaveBSE(
             n_corr_shells=self._n_corr_shells,
             h5_file=os.path.abspath(self._params['bse']['h5_output_file']),
-            bse_info='check',
+            bse_info=bse_info,
             nonlocal_order_parameter=False,
             use_spin_orbit=self._use_spin_orbit,
             beta=self._beta,
@@ -498,9 +506,12 @@ def calc_num_flavors(_ish):
                                                               self._beta, self._n_iw,
                                                               self._sh_quant[ish].Sigma_iw, Gloc_iw_sh[ish],
                                                               self._params['bse']['num_wb'],
-                                                              self._params['bse']['num_wf'], ish, freqs=freqs)
+                                                              self._params['bse']['num_wf'],
+                                                              self._params['bse']['calc_only_chiloc'],
+                                                              ish, freqs=freqs)
 
-            subtract_disconnected(x_loc, g_imp, self.spin_block_names, freqs=freqs)
+            if x_loc is not None:
+                subtract_disconnected(x_loc, g_imp, self.spin_block_names, freqs=freqs)
 
             # Open HDF5 file to improve performance. Close manually.
             bse.h5bse.open('a')
@@ -509,7 +520,8 @@ def calc_num_flavors(_ish):
             for icrsh in range(self._n_corr_shells):
                 if ish == self._sk.corr_to_inequiv[icrsh]:
                     # X_loc
-                    bse.save_xloc(x_loc, icrsh=icrsh)
+                    if x_loc is not None:
+                        bse.save_xloc(x_loc, icrsh=icrsh)
                     # chi_loc
                     if chi_loc is not None:
                         bse.save_chiloc(chi_loc, icrsh=icrsh)
@@ -618,34 +630,39 @@ def dcore_bse(filename, np=1):
     #
     # Construct a parser with default values
     #
-    pars = create_parser()
-
+    pars = create_parser(["model", "system", "impurity_solver", "mpi", "bse"])
     #
     # Parse keywords and store
     #
     pars.read(filename)
-    p = pars.as_dict()
-    parse_parameters(p)
-    seedname = p["model"]["seedname"]
-    p["mpi"]["num_processes"] = np
+    params = pars.as_dict()
+    parse_parameters(params)
+
+    params["mpi"]["num_processes"] = np
+
+    # Delete unnecessary parameters
+    delete_parameters(params, block='model', delete=['interaction', 'density_density', 'kanamori', 'slater_f', 'slater_uj', 'slater_basis', 'interaction_file', 'local_potential_matrix', 'local_potential_factor'])
+    delete_parameters(params, block='model', delete=['bvec'])
+    delete_parameters(params, block='system', retain=['beta', 'n_iw', 'mu', 'fix_mu', 'prec_mu', 'with_dc', 'no_tail_fit'])
+
+    # Summary of input parameters
+    print_parameters(params)
 
     #
     # Load DMFT data
     #
-    solver = DMFTBSESolver(seedname, p, output_file=seedname + '.out.h5')
+    seedname = params["model"]["seedname"]
+    solver = DMFTBSESolver(seedname, params, output_file=seedname + '.out.h5')
     if solver.iteration_number == 0:
         raise RuntimeError("Number of iterations is zero!")
     print("Number of iterations :", solver.iteration_number)
 
     #
-    # Compute data for BSE
+    # Calculate quantities necessary for BSE
     #
     solver.calc_bse()
 
-
-    #
     # Finish
-    #
     print("\n#################  Done  #####################\n")
 
 
diff --git a/src/dcore/dcore_check.py b/src/dcore/dcore_check.py
index af02fcc8..b532c55e 100644
--- a/src/dcore/dcore_check.py
+++ b/src/dcore/dcore_check.py
@@ -95,7 +95,7 @@ def print_chemical_potential(self):
 
     def __plot_init(self):
         if self.plot_called:
-            self.plt.clf()
+            self.plt.close()
             self.plt.figure(figsize=(8, 6))  # default
             return
         self.plot_called = True
diff --git a/src/dcore/dcore_post.py b/src/dcore/dcore_post.py
index 23ef2fb9..db035d89 100644
--- a/src/dcore/dcore_post.py
+++ b/src/dcore/dcore_post.py
@@ -477,7 +477,7 @@ def dcore_post(filename, np=1, prefix=None):
     #
     # Delete unnecessary parameters
     #
-    delete_parameters(p, block='model', delete=['interaction', 'density_density', 'kanamori', 'slater_f', 'slater_uj', 'slater_basis', 'local_potential_matrix', 'local_potential_factor'])
+    delete_parameters(p, block='model', delete=['interaction', 'density_density', 'kanamori', 'slater_f', 'slater_uj', 'slater_basis', 'interaction_file', 'local_potential_matrix', 'local_potential_factor'])
 
     # Summary of input parameters
     print_parameters(p)
diff --git a/src/dcore/dcore_pre.py b/src/dcore/dcore_pre.py
index 57b80a11..b9db75b3 100644
--- a/src/dcore/dcore_pre.py
+++ b/src/dcore/dcore_pre.py
@@ -44,7 +44,12 @@ def __generate_local_potential(p):
     spin_orbit = p["model"]["spin_orbit"]
 
     # read parameters from DFT data
-    skc = SumkDFTCompat(p["model"]["seedname"] + '.h5')
+    try:
+        skc = SumkDFTCompat(p["model"]["seedname"] + '.h5')
+    except RuntimeError as e:
+        print(f"\nERROR: {e}", file=sys.stderr)
+        print("Generate lattice model by running dcore_pre with skip_lattice=False (default) before generating local potential.", file=sys.stderr)
+        sys.exit(1)
 
     assert skc.n_inequiv_shells == n_inequiv_shells
 
@@ -122,7 +127,7 @@ def dcore_pre(input_filenames):
     #
     # Construct a parser with default values
     #
-    pars = create_parser(['model'])
+    pars = create_parser(['model', 'pre'])
     #
     # Parse keywords and store
     #
@@ -140,44 +145,73 @@ def dcore_pre(input_filenames):
     #
     print_parameters(p)
 
+    print("\n@@@@@@@@@@@@@@@@@@@  Generate Model-HDF5 File  @@@@@@@@@@@@@@@@@@@@")
     #
     # remove HDF5 file if exists
     #
     h5_file = p['model']['seedname'] + '.h5'
-    if p['model']['lattice'] != 'external':
-        if os.path.exists(h5_file):
-            print("\nRemoving the existing model HDF5 file...")
+    if os.path.exists(h5_file):
+        print(f"\nFile '{h5_file}' found")
+        if p['pre']['overwrite']:
+            print("  --> overwritten (overwrite=True)")
+        else:
+            print("  --> replaced (overwrite=False)")
             os.remove(h5_file)
+    else:
+        print(f"\nFile '{h5_file}' is created")
+
 
     #
     # Lattice information
     #   -> create h5_file/dft_input
     #
-    print("\n@@@@@@@@@@@@@@@@@@@  Generate Model-HDF5 File  @@@@@@@@@@@@@@@@@@@@\n")
-    lattice_model = create_lattice_model(p)
-    lattice_model.generate_model_file()
+    print("\nGenerating lattice model including H(k)")
+    print(f"  in {h5_file}/dft_input")
+    if p['pre']['skip_lattice']:
+        print("skip")
+    else:
+        lattice_model = create_lattice_model(p)
+        lattice_model.generate_model_file()
 
     #
     # Interaction
     #   -> create h5_file/DCore/umat
     #
     print("\nGenerating U-matrix")
-    generate_umat(p)
+    print(f"  in {h5_file}/DCore/Umat")
+    if p['pre']['skip_umat']:
+        print("skip")
+    else:
+        generate_umat(p)
 
     #
     # Local potential
     #   -> create h5_file/DCore/local_potential
     #
     print("\nGenerating local potential")
-    __generate_local_potential(p)
+    print(f"  in {h5_file}/DCore/LocalPotential")
+    if p['pre']['skip_local_potential']:
+        print("skip")
+    else:
+        __generate_local_potential(p)
 
     #
     # Check HDF5 file
     #
     print('')
     print('@@@@@@@@@@@@@@@@@@@ Check Model-HDF5 file @@@@@@@@@@@@@@@@@@@@')
-    __check_if_Hk_is_hermite(h5_file)
-    print_local_fields(h5_file)
+
+    print("\nChecking H(k)")
+    if p['pre']['skip_lattice']:
+        print("skip")
+    else:
+        __check_if_Hk_is_hermite(h5_file)
+
+    print("\nLocal Fields")
+    if p['pre']['skip_lattice']:
+        print("skip")
+    else:
+        print_local_fields(h5_file)
 
     #
     # Finish
diff --git a/src/dcore/impurity_solvers/alps_cthyb.py b/src/dcore/impurity_solvers/alps_cthyb.py
index 83fb8472..9679a544 100644
--- a/src/dcore/impurity_solvers/alps_cthyb.py
+++ b/src/dcore/impurity_solvers/alps_cthyb.py
@@ -43,7 +43,7 @@ def to_numpy_array(g, block_names):
     Rearrange spins and orbitals so that up and down spins appear alternatingly.
     If there is a single block, we assume that spin and down spins appear alternatignly.
     If there are two blocks, we assume that they are spin1 and spin2 sectors.
-    
+
     Note by HS: The comment above may be wrong. The correct one may be
        If there is a single block, we assume that a up-spin block is followed by a down-spin block.
     """
@@ -125,7 +125,7 @@ def solve(self, rot, mpirun_command, params_kw):
 
         self._solve_impl(rot, mpirun_command, None, params_kw)
 
-    def calc_Xloc_ph(self, rot, mpirun_command, num_wf, num_wb, params_kw):
+    def calc_Xloc_ph(self, rot, mpirun_command, num_wf, num_wb, params_kw, only_chiloc):
         raise RuntimeError("calc_Xloc_ph is not implemented!")
 
     def calc_G2loc_ph_sparse(self, rot, mpirun_command, wsample_ph, params_kw):
diff --git a/src/dcore/impurity_solvers/alps_cthyb_seg.py b/src/dcore/impurity_solvers/alps_cthyb_seg.py
index faf91a7b..11614ae0 100644
--- a/src/dcore/impurity_solvers/alps_cthyb_seg.py
+++ b/src/dcore/impurity_solvers/alps_cthyb_seg.py
@@ -413,29 +413,13 @@ def set_blockgf_from_h5(sigma, group):
         # [(s1,o1), (s2,o2), 0]
         self.quant_to_save['nn_equal_time'] = nn_equal_time[:, :, 0]  # copy
 
-    def calc_Xloc_ph(self, rot, mpirun_command, num_wf, num_wb, params_kw):
+    def calc_Xloc_ph(self, rot, mpirun_command, num_wf, num_wb, params_kw, only_chiloc):
         """
-        compute local G2 in p-h channel
-            X_loc = < c_{i1}^+ ; c_{i2} ; c_{i4}^+ ; c_{i3} >
+        Compute local G2 in p-h channel
 
-        Parameters
-        ----------
-        rot
-        mpirun_command
-        num_wf
-        num_wb
-        params_kw
-
-        Returns
-        -------
-        G2_loc : dict
-            key = (i1, i2, i3, i4)
-            val = numpy.ndarray(n_w2b, 2*n_w2f, 2*n_w2f)
-
-        chi_loc : dict (None if not computed)
-            key = (i1, i2, i3, i4)
-            val = numpy.ndarray(n_w2b)
+        For details, see SolverBase.calc_Xloc_ph
         """
+
         if rot is not None:
             # TODO
             raise NotImplementedError
diff --git a/src/dcore/impurity_solvers/base.py b/src/dcore/impurity_solvers/base.py
index 723c525c..a8a1bae2 100644
--- a/src/dcore/impurity_solvers/base.py
+++ b/src/dcore/impurity_solvers/base.py
@@ -147,29 +147,32 @@ def solve(self, rot, mpirun_command, params):
 
         # Set self.Gimp_iw, self.G_tau, self.Sigma_iw
 
-    def calc_Xloc_ph(self, rot, mpirun_command, num_wf, num_wb, params_kw):
+    def calc_Xloc_ph(self, rot, mpirun_command, num_wf, num_wb, params_kw, only_chiloc):
         """
-        Compute Xloc(m, n, n') in p-h channel
-                and chi_loc(m) (optional)
-
+        Compute local G2 in p-h channel
+            X_loc = < c_{i1}^+ ; c_{i2} ; c_{i4}^+ ; c_{i3} >,  and
+            chi_loc = < c_{i1}^+ c_{i2} ; c_{i4}^+ c_{i3} >  (optional)
 
         Parameters
         ----------
-        num_wf:
+        num_wf: int
             Number of non-negative fermionic frequencies
-        num_wb:
+        num_wb: int
             Number of non-negative bosonic frequencies
+        only_chiloc: bool
+            If True, only chi_loc is computed (no Xloc).
 
         The other parameters are the same as for solve().
 
         Returns
         -------
+        X_loc : dict
+            key : (i1, i2, i3, i4)
+            val : numpy.ndarray(n_w2b, 2*n_w2f, 2*n_w2f)
 
-        Xloc(m, n, n') : 3d ndarray of complex type
-            Data for -num_wf <= n, n' < num_wf and m = 0, 1, ..., num_wb-1.
-
-        chi_loc(m) : 1d ndarray of complex type
-            return None if not computed
+        chi_loc : dict (None if not computed)
+            key : (i1, i2, i3, i4)
+            val : numpy.ndarray(n_w2b)
 
         """
         return NotImplementedError()
diff --git a/src/dcore/impurity_solvers/jo_cthyb_seg.py b/src/dcore/impurity_solvers/jo_cthyb_seg.py
index 8e0c9a2a..1cab4760 100644
--- a/src/dcore/impurity_solvers/jo_cthyb_seg.py
+++ b/src/dcore/impurity_solvers/jo_cthyb_seg.py
@@ -357,28 +357,11 @@ def _read(key):
         # [(s1,o1), (s2,o2), 0]
         # self.quant_to_save['nn_equal_time'] = nn_equal_time[:, :, 0]  # copy
 
-    def calc_Xloc_ph(self, rot, mpirun_command, num_wf, num_wb, params_kw):
+    def calc_Xloc_ph(self, rot, mpirun_command, num_wf, num_wb, params_kw, only_chiloc):
         """
-        compute local G2 in p-h channel
-            X_loc = < c_{i1}^+ ; c_{i2} ; c_{i4}^+ ; c_{i3} >
-
-        Parameters
-        ----------
-        rot
-        mpirun_command
-        num_wf
-        num_wb
-        params_kw
-
-        Returns
-        -------
-        G2_loc : dict
-            key = (i1, i2, i3, i4)
-            val = numpy.ndarray(n_w2b, 2*n_w2f, 2*n_w2f)
-
-        chi_loc : dict (None if not computed)
-            key = (i1, i2, i3, i4)
-            val = numpy.ndarray(n_w2b)
+        Compute local G2 in p-h channel
+
+        For details, see SolverBase.calc_Xloc_ph
         """
         raise NotImplementedError
 
diff --git a/src/dcore/impurity_solvers/pomerol.py b/src/dcore/impurity_solvers/pomerol.py
index cceb6916..508be2cd 100644
--- a/src/dcore/impurity_solvers/pomerol.py
+++ b/src/dcore/impurity_solvers/pomerol.py
@@ -293,44 +293,41 @@ def _read_chiloc(self, params_kw):
         dir_suscep = params_kw.get('dir_suscep', './susceptibility')
         return self._read_common(dir_suscep)
 
-    def calc_Xloc_ph(self, rot, mpirun_command, num_wf, num_wb, params_kw):
+    def calc_Xloc_ph(self, rot, mpirun_command, num_wf, num_wb, params_kw, only_chiloc):
         """
-        compute local G2 in p-h channel
-            X_loc = < c_{i1}^+ ; c_{i2} ; c_{i4}^+ ; c_{i3} >
+        Compute local G2 in p-h channel
 
-        Parameters
-        ----------
-        rot
-        mpirun_command
-        num_wf
-        num_wb
-        params_kw
-
-        Returns
-        -------
-        G2_loc : dict
-            key = (i1, i2, i3, i4)
-            val = numpy.ndarray(n_w2b, 2*n_w2f, 2*n_w2f)
-
-        chi_loc : dict (None if not computed)
-            key = (i1, i2, i3, i4)
-            val = numpy.ndarray(n_w2b)
+        For details, see SolverBase.calc_Xloc_ph
         """
 
-        params_kw['flag_vx'] = 1
-        params_kw['n_w2f'] = num_wf
-        params_kw['n_w2b'] = num_wb
+        # Set parameters
+        if only_chiloc:
+            print("\n Calc only chi_loc (X_loc is not computed)\n")
+            params_kw['flag_vx'] = 0
+        else:
+            params_kw['flag_vx'] = 1
+            params_kw['n_w2f'] = num_wf
+            params_kw['n_w2b'] = num_wb
+
         params_kw['flag_suscep'] = 1
         params_kw['n_wb'] = num_wb
 
+        # Run the impurity solver
         self.solve(rot, mpirun_command, params_kw)
 
-        g2_loc = self._read_g2loc(params_kw)
-        # 1d array --> (wb, wf1, wf2)
-        for key, data in list(g2_loc.items()):
-            g2_loc[key] = data.reshape((num_wb, 2*num_wf, 2*num_wf))
+        # Get results if computed
+        g2_loc = chi_loc = None
+
+        # X_loc
+        if params_kw['flag_vx']:
+            g2_loc = self._read_g2loc(params_kw)
+            # 1d array --> (wb, wf1, wf2)
+            for key, data in list(g2_loc.items()):
+                g2_loc[key] = data.reshape((num_wb, 2*num_wf, 2*num_wf))
 
-        chi_loc = self._read_chiloc(params_kw)
+        # chi_loc
+        if params_kw['flag_suscep']:
+            chi_loc = self._read_chiloc(params_kw)
 
         return g2_loc, chi_loc
 
diff --git a/src/dcore/program_options.py b/src/dcore/program_options.py
index e7114b65..e0ecb78b 100644
--- a/src/dcore/program_options.py
+++ b/src/dcore/program_options.py
@@ -46,7 +46,7 @@ def create_parser(target_sections=None):
     Create a parser for all program options of DCore
     """
     if target_sections is None:
-        parser = TypedParser(['mpi', 'model', 'system', 'impurity_solver', 'control', 'tool', 'post', 'bse', 'vertex', 'sparse_bse'])
+        parser = TypedParser(['mpi', 'model', 'pre', 'system', 'impurity_solver', 'control', 'tool', 'post', 'bse', 'vertex', 'sparse_bse'])
     else:
         parser = TypedParser(target_sections)
 
@@ -84,6 +84,12 @@ def create_parser(target_sections=None):
     parser.add_option("model", "local_potential_matrix", str, "None", "dict of {ish: 'filename'} to specify local potential matrix of ish-th shell")
     parser.add_option("model", "local_potential_factor", str, "1.0", "Prefactors to the local potential matrix (float or list with len=ncor)")
 
+    # [pre]
+    parser.add_option("pre", "skip_lattice", bool, False, "Skip generation of lattice model including H(k) in dcore_pre. If True, data in seedname.h5/dft_input are kept.")
+    parser.add_option("pre", "skip_umat", bool, False, "Skip generation of U matrix in dcore_pre. If True, data in seedname.h5/DCore/Umat are kept.")
+    parser.add_option("pre", "skip_local_potential", bool, False, "Skip generation of local potential in dcore_pre. If True, data in seedname.h5/DCore/LocalPotential are kept.")
+    parser.add_option("pre", "overwrite", bool, True, "Whether seedname.h5 is overwritten. If False, the existing seedname.h5 is repalced.")
+
     # [system]
     parser.add_option("system", "beta", float, 1.0, "Inverse temperature. This parameter is overwritten, if T is given.")
     parser.add_option("system", "T", float, -1.0, "Temperature. If this parameter is given, beta is overwritten by 1/T.")
@@ -179,8 +185,10 @@ def create_parser(target_sections=None):
     parser.add_option("bse", "num_wb", int, 1, "Number of bosonic frequencies (>0)")
     parser.add_option("bse", "num_wf", int, 10, "Number of fermionic frequencies (>0)")
     parser.add_option("bse", "h5_output_file", str, 'dmft_bse.h5', "Output HDF5 file for bse data")
-    parser.add_option("bse", "skip_X0q_if_exists", bool, False, "Skip X_0(q) calc if file already exists")
+    parser.add_option("bse", "skip_X0q_if_exists", bool, False, "[NOT USED] Skip X_0(q) calc if file already exists", OptionStatus.RETIRED)
+    parser.add_option("bse", "skip_X0q", bool, False, "Skip X_0(q) calc")
     parser.add_option("bse", "skip_Xloc", bool, False, "Skip X_loc calc (for RPA)")
+    parser.add_option("bse", "calc_only_chiloc", bool, False, "Calculate only chi_loc but no X_loc (for SCL, rRPA). Do not activate skip_Xloc when using this option.")
     parser.add_option("bse", "use_temp_file", bool, False, "Whether or not temporary file is used in computing X0_q. This option will reduce the memory footprints.")
     parser.add_option("bse", "X0q_qpoints_saved", str, 'quadrant', "Specifies for which q points X0q are saved in a HDF file. quadrant or path to a q_path.dat file.")
 
@@ -297,6 +305,9 @@ def parse_parameters(params):
         if params['post.anacont.spm']['max_iters_opt'] <= 0:
             sys.exit(f"ERROR: max_iters_opt={params['post.anacont.spm']['max_iters_opt']} must be a positive integer.")
         
+    if 'bse' in params:
+        two_options_incompatible(params, ('bse', 'skip_Xloc'), ('bse', 'calc_only_chiloc'))
+
 
 def parse_knode(knode_string):
     """
diff --git a/src/dcore/sumkdft_compat.py b/src/dcore/sumkdft_compat.py
index 936d078e..2cd4689d 100644
--- a/src/dcore/sumkdft_compat.py
+++ b/src/dcore/sumkdft_compat.py
@@ -29,7 +29,7 @@ def read_dft_input_data(file, subgrp, things_to_read):
     values = {}
     with HDFArchive(file, 'r') as ar:
         if not subgrp in ar:
-            raise RuntimeError("subrp " + subgrp + "does not exist in " + file + "!")
+            raise RuntimeError("subgrp " + subgrp + " does not exist in " + file + "!")
         # first read the necessary things:
         for it in things_to_read:
             values[it] = ar[subgrp][it]
diff --git a/src/dcore/sumkdft_opt.py b/src/dcore/sumkdft_opt.py
index 64304f18..d02902e8 100644
--- a/src/dcore/sumkdft_opt.py
+++ b/src/dcore/sumkdft_opt.py
@@ -854,6 +854,7 @@ def eff_atomic_levels(self):
             eff_atlevels = None
         return mpi.bcast(eff_atlevels)
 
+    # This replacement does not work when np_mpi > nk
     def calculate_min_max_band_energies(self):
         # hop = self.hopping
         # diag_hop = numpy.zeros(hop.shape[:-1])
@@ -871,6 +872,16 @@ def calculate_min_max_band_energies(self):
         self.max_band_energy = max_band_energy
         return min_band_energy, max_band_energy
 
+    # Simply set 0
+    def calculate_min_max_band_energies(self):
+        if mpi.is_master_node():
+            warn("Set min_band_energy=0 and max_band_energy=0 when hopping_part is used.")
+        min_band_energy = 0
+        max_band_energy = 0
+        self.min_band_energy = min_band_energy
+        self.max_band_energy = max_band_energy
+        return min_band_energy, max_band_energy
+
     # This method is not used for the moment.
     # Actually, replacement is simple.
     def calc_density_correction(self, filename=None, dm_type='wien2k'):