Merge branch 'develop' into dcore_anacont

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

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")
 
 

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

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)

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

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):