diff --git a/src/dcore/dcore_bse.py b/src/dcore/dcore_bse.py
index dd76096a..efe53609 100644
--- a/src/dcore/dcore_bse.py
+++ b/src/dcore/dcore_bse.py
@@ -49,7 +49,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, only_chiloc, ish, freqs=None):
+                              Sigma_iw, Gloc_iw, num_wb, num_wf, save_chiloc, only_chiloc, ish, freqs=None):
     """
 
     Calculate G2 in an impurity model
@@ -71,6 +71,9 @@ def calc_g2_in_impurity_model(solver_name, solver_params, mpirun_command, basis_
     s_params = copy.deepcopy(solver_params)
     s_params['random_seed_offset'] = 1000 * ish
 
+    s_params['save_chiloc'] = save_chiloc
+    s_params['only_chiloc'] = only_chiloc
+
     work_dir_org = os.getcwd()
     work_dir = 'work/imp_shell' + str(ish) + '_bse'
     if not os.path.isdir(work_dir):
@@ -83,7 +86,7 @@ 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, only_chiloc)
+        xloc, chiloc = sol.calc_Xloc_ph(rot, mpirun_command, num_wf, num_wb, s_params)
     else:
         xloc, chiloc = sol.calc_Xloc_ph_sparse(rot, mpirun_command, freqs, num_wb, s_params)
 
@@ -167,6 +170,56 @@ def g(_i, _j, _w):
                     data[j] -= g(i2, i1, wf1) * g(i3, i4, wf2)
 
 
+def calc_X0_loc(gimp, spin_names, num_wb, num_wf, freqs=None):
+    """
+    calculate X0_loc from G_imp
+        X0_loc[(i1, i2, i3, i4)][wb, wf1, wf2] = - G[(i3, i1)][wf1] * G[(i2, i4)][wf2+wb]
+
+    if freqs is None: box frequency sampling
+    else: sparse sampling
+
+    """
+
+    # g_ij[(i, j)] = data[iw]
+    g_ij = {}
+    for isp, sp in enumerate(spin_names):
+        # data[iw, orb1, orb2]
+        norb = gimp[sp].data.shape[1]
+        assert norb == gimp[sp].data.shape[2]
+        for o1, o2 in product(list(range(norb)), repeat=2):
+            i = o1 + isp*norb
+            j = o2 + isp*norb
+            g_ij[(i, j)] = gimp[sp].data[:, o1, o2]
+
+    assert g_ij[(0, 0)].shape[0] % 2 == 0
+    w0 = g_ij[(0, 0)].shape[0] // 2
+
+    # reduce Matsubara frequencies
+    for key, g in g_ij.items():
+        g_ij[key] = g[w0-num_wf:w0+num_wf+num_wb]
+
+    # get maximum index
+    max_ij = numpy.max(numpy.array([key for key in g_ij.keys()]))
+
+    chi0_loc = {}
+    for i1, i2, i3, i4 in product(range(max_ij+1), repeat=4):
+        # skip if g_ij has no data
+        if (not (i3, i1) in g_ij) or (not (i2, i4) in g_ij):
+            continue
+
+        if freqs is None:
+            # box sampling
+            chi0 = numpy.zeros((num_wb, 2*num_wf), dtype=numpy.complex128)
+            for wb in range(num_wb):
+                chi0[wb, :] = - g_ij[(i3, i1)][0:2*num_wf] * g_ij[(i2, i4)][wb:wb+2*num_wf]
+            chi0_loc[(i1, i2, i3, i4)] = chi0
+        else:
+            # sparse sampling
+            raise NotImplementedError
+
+    return chi0_loc
+
+
 def gen_sparse_freqs_fix_boson(boson_freq, Lambda, sv_cutoff):
     """
     Generate sparse frequency points (wf1, wf2) for fixed bosonic freq wb
@@ -338,6 +391,9 @@ def decompose_index(index):
     def save_xloc(self, xloc_ijkl, icrsh):
         self._save_common(xloc_ijkl, icrsh, 'X_loc')
 
+    def save_x0loc(self, x0loc_ijkl, icrsh):
+        self._save_common(x0loc_ijkl, icrsh, 'X0_loc')
+
     def save_chiloc(self, chiloc_ijkl, icrsh):
         self._save_common(chiloc_ijkl, icrsh, 'chi_loc_in')
 
@@ -424,6 +480,7 @@ def _calc_bse_x0q(self):
         params['div'] = lattice_model.nkdiv()
         params['bse_h5_out_file'] = os.path.abspath(self._params['bse']['h5_output_file'])
         params['use_temp_file'] = self._params['bse']['use_temp_file']
+        params['flag_save_X0_loc'] = self._params['bse']['choice_X0loc'] == 'qsum'
         if self._params['bse']['X0q_qpoints_saved'] == 'quadrant':
             params['X0q_qpoints_saved'] = 'quadrant'
         else:
@@ -507,12 +564,25 @@ def calc_num_flavors(_ish):
                                                               self._sh_quant[ish].Sigma_iw, Gloc_iw_sh[ish],
                                                               self._params['bse']['num_wb'],
                                                               self._params['bse']['num_wf'],
+                                                              self._params['bse']['save_chiloc'],
                                                               self._params['bse']['calc_only_chiloc'],
                                                               ish, freqs=freqs)
 
             if x_loc is not None:
                 subtract_disconnected(x_loc, g_imp, self.spin_block_names, freqs=freqs)
 
+            if self._params['bse']['choice_X0loc'] == 'imp':
+                print("calc_X0_loc start")
+                params = dict(
+                    spin_names = self.spin_block_names,
+                    num_wf = self._params['bse']['num_wf'],
+                    num_wb = self._params['bse']['num_wb'],
+                )
+                x0_loc = calc_X0_loc(g_imp, **params, freqs=freqs)
+                print("calc_X0_loc end")
+            else:
+                x0_loc = None
+
             # Open HDF5 file to improve performance. Close manually.
             bse.h5bse.open('a')
 
@@ -525,6 +595,9 @@ def calc_num_flavors(_ish):
                     # chi_loc
                     if chi_loc is not None:
                         bse.save_chiloc(chi_loc, icrsh=icrsh)
+                    # X0_loc
+                    if x0_loc is not None:
+                        bse.save_x0loc(x0_loc, icrsh=icrsh)
 
             bse.h5bse.close()
 
@@ -532,8 +605,12 @@ def calc_bse(self):
         """
         Compute data for BSE
         """
-        self._calc_bse_x0q()
-        self._calc_bse_xloc()
+        self._calc_bse_x0q()  # X_0(q)
+        self._calc_bse_xloc()  # X_{loc}
+
+        # X_{0,loc} is computed in
+        #   _calc_bse_x0q()  if calc_X0loc_from_X0q = True (default)
+        #   _calc_bse_xloc() if                     = False
 
 
     def fit_Xloc(self, save_only):
diff --git a/src/dcore/impurity_solvers/alps_cthyb.py b/src/dcore/impurity_solvers/alps_cthyb.py
index 9679a544..a43ea0a5 100644
--- a/src/dcore/impurity_solvers/alps_cthyb.py
+++ b/src/dcore/impurity_solvers/alps_cthyb.py
@@ -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, only_chiloc):
+    def calc_Xloc_ph(self, rot, mpirun_command, num_wf, num_wb, params_kw):
         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 11614ae0..6c20c242 100644
--- a/src/dcore/impurity_solvers/alps_cthyb_seg.py
+++ b/src/dcore/impurity_solvers/alps_cthyb_seg.py
@@ -23,6 +23,7 @@
 from dcore._dispatcher import *
 from ..tools import make_block_gf, launch_mpi_subprocesses, extract_H0, umat2dd, get_block_size, expand_path
 from .base import SolverBase, rotate_basis
+from dcore.program_options import parse_save_chiloc
 
 
 def to_numpy_array(g, names):
@@ -413,23 +414,20 @@ 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, only_chiloc):
+    def calc_Xloc_ph(self, rot, mpirun_command, num_wf, num_wb, params_kw):
         """
         Compute local G2 in p-h channel
 
         For details, see SolverBase.calc_Xloc_ph
         """
 
-        if rot is not None:
-            # TODO
-            raise NotImplementedError
-
-        use_chi_loc = False
+        save_chiloc = parse_save_chiloc(params_kw['save_chiloc'], default=False)
+        only_chiloc = params_kw['only_chiloc']
 
         params_kw['cthyb.MEASURE_g2w'] = 1
         params_kw['cthyb.N_w2'] = num_wf
         params_kw['cthyb.N_W'] = num_wb
-        if use_chi_loc:
+        if save_chiloc:
             params_kw['cthyb.MEASURE_nnw'] = 1
 
         self.solve(rot, mpirun_command, params_kw)
@@ -476,7 +474,7 @@ def get_g2(_i, _j, _wb, _wf1, _wf2):
         # Save chi(wb)
         # [(s1,o1), (s2,o2), wb]
         chi_dict = None
-        if use_chi_loc:
+        if save_chiloc:
             chi_re = self._get_results("nnw_re", num_wb, orbital_symmetrize=True)
             chi_im = self._get_results("nnw_im", num_wb, orbital_symmetrize=True)
             chi_loc = chi_re + chi_im * 1.0J
@@ -487,6 +485,13 @@ def get_g2(_i, _j, _wb, _wf1, _wf2):
             for i1, i2 in product(range(2*self.n_orb), repeat=2):
                 chi_dict[(i1, i1, i2, i2)] = chi_loc[i1, i2]
 
+        # Rotate g2_dict and chi_dict back to the original basis
+        if rot is not None:
+            print("Error: 'basis_rotation' for two-particle quantities is currently not supported.", file=sys.stderr)
+            sys.exit(1)
+
+            rotate_basis(rot, self.use_spin_orbit, None, direction='backward', X_dict=g2_dict, chi_dict=chi_dict)
+
         return g2_dict, chi_dict
 
     def calc_G2loc_ph_sparse(self, rot, mpirun_command, freqs_ph, num_wb, params_kw):
diff --git a/src/dcore/impurity_solvers/base.py b/src/dcore/impurity_solvers/base.py
index a8a1bae2..0ebe7ece 100644
--- a/src/dcore/impurity_solvers/base.py
+++ b/src/dcore/impurity_solvers/base.py
@@ -147,7 +147,7 @@ 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, only_chiloc):
+    def calc_Xloc_ph(self, rot, mpirun_command, num_wf, num_wb, params_kw):
         """
         Compute local G2 in p-h channel
             X_loc = < c_{i1}^+ ; c_{i2} ; c_{i4}^+ ; c_{i3} >,  and
@@ -159,10 +159,12 @@ def calc_Xloc_ph(self, rot, mpirun_command, num_wf, num_wb, params_kw, only_chil
             Number of non-negative fermionic frequencies
         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().
+        params_kw includes the following parameters in addition to the ones used in solve().
+
+            'only_chiloc': bool, if True, only chi_loc is computed (no Xloc).
+            'save_chiloc': bool, if True, chi_loc is saved.
 
         Returns
         -------
@@ -259,7 +261,7 @@ def is_Floc_computable(cls):
 
 
 
-def rotate_basis(rot, use_spin_orbit, u_matrix, Gfs=[], direction='forward'):
+def rotate_basis(rot, use_spin_orbit, u_matrix, Gfs=[], direction='forward', X_dict=None, chi_dict=None):
     """
     Rotate all Gf-like objects and U-matrix to the basis defined by rot
 
@@ -268,17 +270,17 @@ def rotate_basis(rot, use_spin_orbit, u_matrix, Gfs=[], direction='forward'):
     """
 
     if direction == 'forward':
-        return _rotate_basis(rot, u_matrix, use_spin_orbit, Gfs)
+        return _rotate_basis(rot, u_matrix, use_spin_orbit, Gfs, X_dict, chi_dict)
     elif direction == 'backward':
         rot_conj_trans = {}
         for name, r in list(rot.items()):
             rot_conj_trans[name] = r.conjugate().transpose()
-        return _rotate_basis(rot_conj_trans, u_matrix, use_spin_orbit, Gfs)
+        return _rotate_basis(rot_conj_trans, u_matrix, use_spin_orbit, Gfs, X_dict, chi_dict)
     else:
         raise RuntimeError("Unknown direction " + direction)
 
 
-def _rotate_basis(rot, u_matrix, use_spin_orbit, Gfs):
+def _rotate_basis(rot, u_matrix, use_spin_orbit, Gfs, X_dict, chi_dict):
     """
     Rotate all Gf-like object and U matrix to a new local basis defined by "rot".
     :param rot: matrix
@@ -298,11 +300,52 @@ def _rotate_basis(rot, u_matrix, use_spin_orbit, Gfs):
         for bname, gf in G:
             gf.from_L_G_R(rot[bname].transpose().conjugate(), gf, rot[bname])
 
+    n_flavors = rot_spin_full.shape[0]
+
+    if X_dict is not None:
+        shape = next(iter(X_dict.values())).shape  # (num_wb, num_wf, num_wf)
+        assert len(shape) == 3
+        array = numpy.zeros((n_flavors, n_flavors, n_flavors, n_flavors) + shape, dtype=numpy.complex128)
+        _set_from_dict(X_dict, array)
+
+        # X_{1234}(iW, iw, iw') = < c_1^+(iw) c_2(iw+iW) c_4^+(iw'+iW) c_3(iw') >
+        array = numpy.einsum("ijklxyz,im,jn,ko,lp -> mnopxyz", array,
+                    rot_spin_full, numpy.conj(rot_spin_full), numpy.conj(rot_spin_full), rot_spin_full)
+
+        X_dict.clear()
+        _set_to_dict(X_dict, array)
+
+    if chi_dict is not None:
+        shape = next(iter(chi_dict.values())).shape  # (num_wb,)
+        assert len(shape) == 1
+        array = numpy.zeros((n_flavors, n_flavors, n_flavors, n_flavors) + shape, dtype=numpy.complex128)
+        _set_from_dict(chi_dict, array)
+
+        # chi_{1234}(iW) = < c_1^+ c_2 c_4^+ c_3 >(iW)
+        array = numpy.einsum("ijklx,im,jn,ko,lp -> mnopx", array,
+                    rot_spin_full, numpy.conj(rot_spin_full), numpy.conj(rot_spin_full), rot_spin_full)
+
+        chi_dict.clear()
+        _set_to_dict(chi_dict, array)
+
     if not u_matrix is None:
         return numpy.einsum("ijkl,im,jn,ko,lp", u_matrix,
                                     numpy.conj(rot_spin_full), numpy.conj(rot_spin_full), rot_spin_full, rot_spin_full)
 
 
+def _set_from_dict(x_dict, x_array):
+    for key, val in x_dict.items():
+        i, j, k, l = key
+        x_array[i, j, k, l] = val  # val is np.array
+
+
+def _set_to_dict(x_dict, x_array):
+    n1, n2, n3, n4 = x_array.shape[:4]
+    for i, j, k, l in product(range(n1), range(n2), range(n3), range(n4)):
+        if not numpy.all(np.abs(x_array[i, j, k, l]) < 1e-8):  # if not zero matrix
+            x_dict[(i, j, k, l)] = x_array[i, j, k, l]
+
+
 class PytriqsMPISolver(SolverBase):
 
     def __init__(self, beta, gf_struct, u_mat, n_iw=1025):
diff --git a/src/dcore/impurity_solvers/jo_cthyb_seg.py b/src/dcore/impurity_solvers/jo_cthyb_seg.py
index 1cab4760..ac480fc9 100644
--- a/src/dcore/impurity_solvers/jo_cthyb_seg.py
+++ b/src/dcore/impurity_solvers/jo_cthyb_seg.py
@@ -24,6 +24,7 @@
 from dcore._dispatcher import *
 from ..tools import make_block_gf, launch_mpi_subprocesses, extract_H0, umat2dd, get_block_size, expand_path
 from .base import SolverBase, rotate_basis
+from dcore.program_options import parse_save_chiloc
 
 
 def to_numpy_array(g, names):
@@ -352,86 +353,127 @@ def _read(key):
             rotate_basis(rot, self.use_spin_orbit, None, [self._Sigma_iw, self._Gimp_iw], direction='backward')
             # rotate_basis(rot, self.use_spin_orbit, None, self._Gimp_iw, direction='backward')
 
-        #   self.quant_to_save['nn_equal_time']
-        # nn_equal_time =
-        # [(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, only_chiloc):
+    def calc_Xloc_ph(self, rot, mpirun_command, num_wf, num_wb, params_kw):
         """
         Compute local G2 in p-h channel
 
         For details, see SolverBase.calc_Xloc_ph
         """
-        raise NotImplementedError
 
-        if rot is not None:
-            # TODO
-            raise NotImplementedError
+        save_chiloc = parse_save_chiloc(params_kw['save_chiloc'], default=False)
+        only_chiloc = params_kw['only_chiloc']
 
-        use_chi_loc = False
+        if save_chiloc:
+            print("Warning: [bse]save_chiloc=True. The transverse spin susceptibility and some orbital susceptibilities calculated via BSE will be incorrect, although the longitudinal spin susceptibility shows faster convergence with respect to num_wf. [bse]save_chiloc=False is recommended. Note, however, that the spin rotational symmetry remains broken due to the use of density-density interactions.", file=sys.stderr)
 
-        params_kw['cthyb.MEASURE_g2w'] = 1
-        params_kw['cthyb.N_w2'] = num_wf
-        params_kw['cthyb.N_W'] = num_wb
-        if use_chi_loc:
-            params_kw['cthyb.MEASURE_nnw'] = 1
+        params_kw['control.flag_tp'] = 'true'
+        params_kw['control.n_tp'] = 2**(int(np.log2(self.n_iw)) - 4)  # TODO
+        params_kw['control.n_tp2'] = 2**(int(np.log2(self.n_iw)) - 1)
+
+        if not only_chiloc:
+            params_kw['control.flag_vx'] = 'true'
+            params_kw['control.n_vx1'] = num_wf
+            params_kw['control.n_vx2'] = num_wb
 
         self.solve(rot, mpirun_command, params_kw)
 
-        # Save G2(wb, wf, wf')
-        # [(s1,o1), (s2,o2), (wb,wf,wf')]
-        g2_re = self._get_results("g2w_re", 4*num_wf*num_wf*num_wb, orbital_symmetrize=False)
-        g2_im = self._get_results("g2w_im", 4*num_wf*num_wf*num_wb, orbital_symmetrize=False)
-        g2_loc = (g2_re + g2_im * 1.0J) / self.beta
-        g2_loc = g2_loc.reshape((2*self.n_orb, 2*self.n_orb) + (num_wb, 2*num_wf, 2*num_wf))
-        # assign to dict
-        g2_dict = {}
-        for i1, i2 in product(range(2*self.n_orb), repeat=2):
-            g2_dict[(i1, i1, i2, i2)] = g2_loc[i1, i2]
-
-        # return g2_loc for arbitrary wb including wb<0
-        def get_g2(_i, _j, _wb, _wf1, _wf2):
-            try:
-                if _wb >= 0:
-                    return g2_loc[_i, _j, _wb, _wf1, _wf2]
-                else:
-                    # G2_iijj(wb, wf, wf') = G2_jjii(-wb, -wf', -wf)^*
-                    return numpy.conj(g2_loc[_j, _i, -_wb, -(1+_wf2), -(1+_wf1)])
-            except IndexError:
-                return 0
-
-        # Convert G2_iijj -> G2_ijij
-        g2_loc_tr = numpy.zeros(g2_loc.shape, dtype=complex)
-        for i1, i2 in product(range(2*self.n_orb), repeat=2):
-            for wb in range(num_wb):
-                for wf1, wf2 in product(range(2 * num_wf), repeat=2):
-                    # G2_ijij(wb, wf, wf') = -G2_iijj(wf-wf', wf'+wb, wf')^*
-                    g2_loc_tr[i1, i2, wb, wf1, wf2] = -get_g2(i1, i2, wf1-wf2, wf2+wb, wf2)
-        # assign to dict
-        for i1, i2 in product(range(2*self.n_orb), repeat=2):
-            # exclude i1=i2, which was already assigned by g2_loc
-            if i1 != i2:
-                g2_dict[(i1, i2, i1, i2)] = g2_loc_tr[i1, i2]
-
-        # Occupation number
-        # [(s1,o1)]
-        occup = self._get_occupation()
-
-        # Save chi(wb)
-        # [(s1,o1), (s2,o2), wb]
+        # Get G2
+        if only_chiloc:
+            g2_dict = None
+        else:
+            # Gf[wf, i]
+            # Read Gf
+            data = numpy.loadtxt("wmeasure_Gf_w.dat")
+            gf = data[:, 1::2] + 1j * data[:, 2::2]
+            assert gf.shape == (2*num_wf + num_wb, 2*self.n_orb)
+
+            # x0_lo[i, j, wf, wf']  (wb=0)
+            # x0_{12,34}(0; wf, wf') = gf_{21}(iw) * gf_{34}(iw')
+            # x0_{11,33}(0; wf, wf') = gf_{11}(iw) * gf_{33}(iw')  [Gf is diagonal]
+            x0_lo = np.zeros((2*self.n_orb, 2*self.n_orb, 2*num_wf, 2*num_wf), dtype=numpy.complex128)
+            for i1, i3 in product(range(2*self.n_orb), repeat=2):
+                x0_lo[i1, i3, :, :] = gf[:2*num_wf, None, i1] * gf[None, :2*num_wf, i3]
+
+            # x0_tr[i, j, wb, wf]  (wf=wf')
+            # x0_{12,34}(wb; wf) = - gf_{31}(wf) * gf_{24}(wf+wb)
+            # x0_{12,12}(wb; wf) = - gf_{11}(wf) * gf_{22}(wf+wb)  [Gf is diagonal]
+            x0_tr = np.zeros((2*self.n_orb, 2*self.n_orb, num_wb, 2*num_wf), dtype=numpy.complex128)
+            for i1, i2 in product(range(2*self.n_orb), repeat=2):
+                for wb in range(num_wb):
+                    x0_tr[i1, i2, wb, :] = - gf[0:2*num_wf, i1] * gf[wb:wb+2*num_wf, i2]
+
+            # Read gamma
+            def read_vertex(filename):
+                data = numpy.loadtxt(filename)
+                gamma = data[:, 6::2] + 1j * data[:, 7::2]
+                assert gamma.shape == (num_wb * 2*num_wf * 2*num_wf, 2*self.n_orb * 2*self.n_orb)
+                gamma = gamma.reshape(2*num_wf, 2*num_wf, num_wb, 2*self.n_orb, 2*self.n_orb)
+                # (wf, wf', wb, i, j) -> (i, j, wb, wf, wf')
+                gamma = gamma.transpose([3, 4, 2, 0, 1])
+                # (i, j, wb, wf, wf')
+                assert gamma.shape == (2*self.n_orb, 2*self.n_orb, num_wb, 2*num_wf, 2*num_wf)
+                return -gamma / self.beta
+
+            # gamma[i, j, wb, wf, wf']
+            gamma_lo = read_vertex("vertex_lo.dat")
+            gamma_tr = read_vertex("vertex_tr.dat")
+
+            # gamma -> connected G2
+            # G2[i, j, wb, wf, wf']
+            #  = x0_{12}(wb; wf) * gamma_{12,34}(wb; wf, wf') * x0_{34}(wb; wf')
+
+            # for lo
+            # x0_{11}(wb; wf) * gamma_{11,33}(wb; wf, wf') * x0_{33}(wb; wf')
+            g2_lo = np.einsum("iixy, ijxyz, jjxz -> ijxyz", x0_tr, gamma_lo, x0_tr)
+
+            # for tr
+            # x0_{12}(wb; wf) * gamma_{12,12}(wb; wf, wf') * x0_{12}(wb; wf')
+            g2_tr = np.einsum("ijxy, ijxyz, ijxz -> ijxyz", x0_tr, gamma_tr, x0_tr)
+
+            # Sum disconnected part
+            # for lo
+            # wb == 0
+            g2_lo[:, :, 0, :, :] += x0_lo[:, :, :, :]
+            # i == j, wf == wf'
+            for i in range(2*self.n_orb):
+                for wf in range(2*num_wf):
+                    g2_lo[i, i, :, wf, wf] += x0_tr[i, i, :, wf]
+
+            # for tr (i != j)
+            # wf == wf'
+            for wf in range(2*num_wf):
+                g2_tr[:, :, :, wf, wf] += x0_tr[:, :, :, wf]
+
+            g2_dict = {}
+            for i1, i2 in product(range(2*self.n_orb), repeat=2):
+                g2_dict[(i1, i1, i2, i2)] = g2_lo[i1, i2]
+                if i1 != i2:
+                    g2_dict[(i1, i2, i1, i2)] = g2_tr[i1, i2]
+
+        # Get chi
         chi_dict = None
-        if use_chi_loc:
-            chi_re = self._get_results("nnw_re", num_wb, orbital_symmetrize=True)
-            chi_im = self._get_results("nnw_im", num_wb, orbital_symmetrize=True)
-            chi_loc = chi_re + chi_im * 1.0J
-            # subtract <n><n>
-            chi_loc[:, :, 0] -= occup[:, None] * occup[None, :] * self.beta
-            # assign to dict
+        if save_chiloc:
+            # Read chi_loc
+            data = numpy.loadtxt("chi_w.dat")
+            chi_loc = data[:num_wb, 3::2] + 1j * data[:num_wb, 4::2]
+            assert chi_loc.shape == (num_wb, (2*self.n_orb)**2)
+            chi_loc = chi_loc.reshape((num_wb, 2*self.n_orb, 2*self.n_orb))
+            # (wb, i, j) -> (i, j, wb)
+            chi_loc = chi_loc.transpose([1, 2, 0])
+            assert chi_loc.shape == (2*self.n_orb, 2*self.n_orb, num_wb)
+
             chi_dict = {}
             for i1, i2 in product(range(2*self.n_orb), repeat=2):
                 chi_dict[(i1, i1, i2, i2)] = chi_loc[i1, i2]
 
+        # Rotate g2_dict and chi_dict back to the original basis
+        if rot is not None:
+            print("Error: 'basis_rotation' for two-particle quantities is currently not supported.", file=sys.stderr)
+            sys.exit(1)
+
+            # TODO: rotate_basis has been implemented, but further test and debugging are necessary.
+            rotate_basis(rot, self.use_spin_orbit, None, direction='backward', X_dict=g2_dict, chi_dict=chi_dict)
+
         return g2_dict, chi_dict
 
     def calc_G2loc_ph_sparse(self, rot, mpirun_command, freqs_ph, num_wb, params_kw):
diff --git a/src/dcore/impurity_solvers/pomerol.py b/src/dcore/impurity_solvers/pomerol.py
index 4de4067a..a0464391 100644
--- a/src/dcore/impurity_solvers/pomerol.py
+++ b/src/dcore/impurity_solvers/pomerol.py
@@ -25,6 +25,7 @@
 
 from ..tools import make_block_gf, launch_mpi_subprocesses, extract_H0, extract_bath_params, expand_path
 from .base import SolverBase
+from dcore.program_options import parse_save_chiloc
 
 VERSION_REQUIRED = 1.5
 
@@ -311,13 +312,16 @@ 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, only_chiloc):
+    def calc_Xloc_ph(self, rot, mpirun_command, num_wf, num_wb, params_kw):
         """
         Compute local G2 in p-h channel
 
         For details, see SolverBase.calc_Xloc_ph
         """
 
+        save_chiloc = parse_save_chiloc(params_kw['save_chiloc'], default=True)
+        only_chiloc = params_kw['only_chiloc']
+
         # Set parameters
         if only_chiloc:
             print("\n Calc only chi_loc (X_loc is not computed)\n")
@@ -327,8 +331,9 @@ def calc_Xloc_ph(self, rot, mpirun_command, num_wf, num_wb, params_kw, only_chil
             params_kw['n_w2f'] = num_wf
             params_kw['n_w2b'] = num_wb
 
-        params_kw['flag_suscep'] = 1
-        params_kw['n_wb'] = num_wb
+        if save_chiloc:
+            params_kw['flag_suscep'] = 1
+            params_kw['n_wb'] = num_wb
 
         # Run the impurity solver
         self.solve(rot, mpirun_command, params_kw)
@@ -344,7 +349,8 @@ def calc_Xloc_ph(self, rot, mpirun_command, num_wf, num_wb, params_kw, only_chil
                 g2_loc[key] = data.reshape((num_wb, 2*num_wf, 2*num_wf))
 
         # chi_loc
-        if params_kw['flag_suscep']:
+        # if params_kw['flag_suscep']:
+        if save_chiloc:
             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 b780f42c..993fb079 100644
--- a/src/dcore/program_options.py
+++ b/src/dcore/program_options.py
@@ -48,7 +48,7 @@ def create_parser(target_sections=None):
         parser = TypedParser(['mpi', 'model', 'pre', 'system', 'impurity_solver', 'control', 'post', 'post.anacont', 'post.anacont.pade', 'post.anacont.spm', 'post.anacont.spm.solver', 'post.spectrum', 'post.check', 'bse', 'vertex', 'sparse_bse'])
     else:
         parser = TypedParser(target_sections)
-    
+
     # tool is removed but read for warning
     parser.allow_undefined_options("tool")
     parser.allow_undefined_options("post.anacont.spm.solver")
@@ -180,9 +180,11 @@ def create_parser(target_sections=None):
     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", "save_chiloc", str, 'default', "Chosen from 'default', 'True', or 'False'. If True, chi_loc is saved to get better convergence against num_wf in ChiQ. But False is recommended for CTHYB segment algorithm. For 'default', True is set for pomerol solver, while False is set for JO/cthyb-seg and ALPS/cthyb-seg solvers.")
     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.")
+    parser.add_option("bse", "choice_X0loc", str, 'qsum', "How to compute X0_loc. 'qsum': by taking q-sum of X_0(q); 'imp': from G_imp. These are not equivalent if G_loc and G_imp are not equal.")
 
     return parser
 
@@ -203,6 +205,12 @@ def two_options_incompatible(params, option1, option2):
         sys.exit(f"ERROR: [{block1}]{param1} and [{block2}]{param2} cannot be True at the same time.")
 
 
+def check_option_choices(params, option, choices):
+    block, param = option
+    if params[block][param] not in choices:
+        sys.exit(f"ERROR: [{block}]{param} must be chosen from {choices}. {repr(params[block][param])} was given.")
+
+
 def parse_parameters(params):
     """
     Parse some parameters in a parameter
@@ -329,6 +337,9 @@ def parse_parameters(params):
     if 'bse' in params:
         two_options_incompatible(params, ('bse', 'skip_Xloc'), ('bse', 'calc_only_chiloc'))
 
+        check_option_choices(params, ('bse', 'save_chiloc'), ['default', 'True', 'False'])
+        check_option_choices(params, ('bse', 'choice_X0loc'), ['qsum', 'imp'])
+
 
 def parse_knode(knode_string):
     """
@@ -389,6 +400,15 @@ def parse_bvec(bvec_string):
     return bvec
 
 
+def parse_save_chiloc(save_chiloc_in, default : bool):
+    assert save_chiloc_in in ('default', 'True', 'False')
+    if save_chiloc_in == 'default':
+        save_chiloc = default
+    else:
+        save_chiloc = save_chiloc_in == 'True'  # convert str to bool
+    return save_chiloc
+
+
 def delete_parameters(p, block, delete=None, retain=None):
     """
     Delete parameters
diff --git a/src/dcore/sumkdft_workers/bse_worker.py b/src/dcore/sumkdft_workers/bse_worker.py
index 42f36c78..c9dd1907 100644
--- a/src/dcore/sumkdft_workers/bse_worker.py
+++ b/src/dcore/sumkdft_workers/bse_worker.py
@@ -35,6 +35,7 @@ def run(self):
 
         sk.save_X0q_for_bse(list_wb=self.params['list_wb'],
                             n_wf_cutoff=self.params['n_wf_G2'],
+                            flag_save_X0_loc=self.params['flag_save_X0_loc'],
                             qpoints_saved=self.params['X0q_qpoints_saved'],
                             h5_file=self.params['bse_h5_out_file'],
                             temp_file=temp_file,