Merge pull request #142 from issp-center-dev/bse_only_chiloc

[bse]calc_only_chiloc option

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, delete_parameters, print_parameters
 from dcore.tools import *
 from dcore import impurity_solvers
 from .sumkdft_workers.launcher import run_sumkdft
@@ -52,7 +52,7 @@ def compare_str_list(list1, list2):
 
 
 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 +86,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)
@@ -498,9 +503,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 +517,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,7 +627,7 @@ 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
@@ -629,6 +638,13 @@ def dcore_bse(filename, np=1):
     seedname = p["model"]["seedname"]
     p["mpi"]["num_processes"] = np
 
+    # Delete unnecessary parameters
+    delete_parameters(p, block='model', delete=['bvec'])
+    delete_parameters(p, block='system', retain=['beta', 'n_iw', 'mu', 'fix_mu', 'prec_mu', 'with_dc', 'no_tail_fit'])
+
+    # Summary of input parameters
+    print_parameters(p)
+
     #
     # Load DMFT data
     #
@@ -638,14 +654,11 @@ def dcore_bse(filename, np=1):
     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/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):

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

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

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
 

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
 

src/dcore/program_options.py

@@ -150,6 +150,7 @@ def create_parser(target_sections=None):
     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_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).")
     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.")
 
@@ -240,6 +241,9 @@ def parse_parameters(params):
         if params['tool']['n_pade_max'] < 0:
             params['tool']['n_pade_max'] = params['system']['n_iw']
 
+    if 'bse' in params:
+        two_options_incompatible(params, ('bse', 'skip_Xloc'), ('bse', 'calc_only_chiloc'))
+
 
 def parse_knode(knode_string):
     """

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