diff --git a/src/hwave/solver/rpa.py b/src/hwave/solver/rpa.py
index c6aff92..5657cec 100644
--- a/src/hwave/solver/rpa.py
+++ b/src/hwave/solver/rpa.py
@@ -51,30 +51,15 @@ class Lattice:
     Note:
       Bx, By, Bz must be chosen so that Lx, Ly, Lz are divisable 
       by Bx, By, Bz, respectively. otherwise, initialization fails.
-    """
 
+    """
     def __init__(self, params):
-        """
-        Initialize the Lattice object.
-
-        Args:
-            params (dict): Dictionary containing lattice parameters.
-        """
         logger.debug(">>> Lattice.__init__")
 
         self._init_lattice(params)
         self._show_params()
 
     def _init_lattice(self, params):
-        """
-        Initialize lattice parameters from the given dictionary.
-
-        Args:
-            params (dict): Dictionary containing lattice parameters.
-
-        Raises:
-            SystemExit: If required parameters are missing or invalid.
-        """
         logger.debug(">>> Lattice._init_lattice")
 
         if "CellShape" not in params:
@@ -133,9 +118,6 @@ def _init_lattice(self, params):
         self.nvol = nvol
 
     def _show_params(self):
-        """
-        Log the lattice parameters.
-        """
         logger.info("Lattice parameters:")
         logger.info("    CellShape       = {}".format(self.cellshape))
         logger.info("    cell volume     = {}".format(self.cellvol))
@@ -151,14 +133,6 @@ class Interaction:
     Construct Hamiltonian from input
     """
     def __init__(self, lattice, param_ham, info_mode):
-        """
-        Initialize the Interaction object.
-
-        Args:
-            lattice (Lattice): The lattice object.
-            param_ham (dict): Dictionary containing Hamiltonian parameters.
-            info_mode (dict): Dictionary containing mode information.
-        """
         logger.debug(">>> Interaction.__init__")
 
         self.lattice = lattice
@@ -169,6 +143,7 @@ def __init__(self, lattice, param_ham, info_mode):
         self._has_interaction_exchange = False
         self._has_interaction_pairhop = False
 
+
         # mode options
         self.enable_spin_orbital = info_mode.get("enable_spin_orbital", False)
 
@@ -184,45 +159,18 @@ def __init__(self, lattice, param_ham, info_mode):
         pass
 
     def has_interaction(self):
-        """
-        Check if there is any interaction.
-
-        Returns:
-            bool: True if there is any interaction, False otherwise.
-        """
         return self._has_interaction
 
     def has_interaction_coulomb(self):
-        """
-        Check if there is Coulomb interaction.
-
-        Returns:
-            bool: True if there is Coulomb interaction, False otherwise.
-        """
         return self._has_interaction_coulomb
 
     def has_interaction_exchange(self):
-        """
-        Check if there is exchange interaction.
-
-        Returns:
-            bool: True if there is exchange interaction, False otherwise.
-        """
         return self._has_interaction_exchange
 
     def has_interaction_pairhop(self):
-        """
-        Check if there is pair hopping interaction.
-
-        Returns:
-            bool: True if there is pair hopping interaction, False otherwise.
-        """
         return self._has_interaction_pairhop
 
     def _init_interaction(self):
-        """
-        Initialize interaction parameters.
-        """
         logger.debug(">>> Interaction._init_interaction")
 
         # reinterpret interaction coefficient on sublattice
@@ -246,15 +194,6 @@ def _init_interaction(self):
         pass
 
     def _reshape_geometry(self, geom):
-        """
-        Reshape the geometry parameters for the sublattice.
-
-        Args:
-            geom (dict): Dictionary containing geometry parameters.
-
-        Returns:
-            dict: Reshaped geometry parameters.
-        """
         logger.debug(">>> Interaction._reshape_geometry")
 
         Bx,By,Bz = self.lattice.subshape
@@ -280,16 +219,6 @@ def _reshape_geometry(self, geom):
         return geom_new
 
     def _reshape_interaction(self, ham, enable_spin_orbital):
-        """
-        Reshape the interaction parameters for the sublattice.
-
-        Args:
-            ham (dict): Dictionary containing Hamiltonian parameters.
-            enable_spin_orbital (bool): Flag to enable spin-orbital interaction.
-
-        Returns:
-            dict: Reshaped Hamiltonian parameters.
-        """
         logger.debug(">>> Interaction._reshape_interaction")
 
         Bx,By,Bz = self.lattice.subshape
@@ -319,7 +248,7 @@ def _round(x, n):
 
             for bz,by,bx in itertools.product(range(Bz),range(By),range(Bx)):
 
-                # original cell index of both ends
+                # original cell index of both endes
                 #   x0 -> x1=x0+r
                 x0,y0,z0 = bx, by, bz
                 x1,y1,z1 = x0 + rx, y0 + ry, z0 + rz
@@ -344,13 +273,6 @@ def _round(x, n):
         return ham_new
 
     def _export_interaction(self, type, file_name):
-        """
-        Export the interaction parameters to a file.
-
-        Args:
-            type (str): Type of interaction.
-            file_name (str): Name of the file to save the interaction parameters.
-        """
         logger.debug(">>> Interaction._export_interaction")
 
         intr = self.param_ham[type]
@@ -385,16 +307,6 @@ def _export_interaction(self, type, file_name):
                     ))
 
     def _make_ham_trans(self):
-        """
-        Create the Hamiltonian transfer terms and perform Fourier transform.
-
-        This method initializes the Hamiltonian transfer terms from the parameters,
-        reshapes them, and performs a Fourier transform to obtain the Hamiltonian
-        in momentum space.
-
-        Raises:
-            Warning: If 'Transfer' key is not found in the Hamiltonian parameters.
-        """
         logger.debug(">>> Interaction._make_ham_trans")
 
         nx,ny,nz = self.lattice.shape
@@ -411,7 +323,7 @@ def _make_ham_trans(self):
             self.ham_extern_q = None
             return
 
-        if self.enable_spin_orbital:
+        if self.enable_spin_orbital == True:
             # assume orbital index includes spin index
             tab_r = np.zeros((nx,ny,nz,nd,nd), dtype=np.complex128)
 
@@ -471,20 +383,9 @@ def _make_ham_trans(self):
         else:
             self.ham_extern_r = None
             self.ham_extern_q = None
-    def _make_ham_inter(self):
-        """
-        Create the interaction Hamiltonian terms and perform Fourier transform.
 
-        This method initializes the interaction Hamiltonian terms from the parameters,
-        reshapes them, and performs a Fourier transform to obtain the Hamiltonian
-        in momentum space.
 
-        The interaction Hamiltonian includes Coulomb, Hund, Ising, PairLift, Exchange,
-        and PairHop interactions.
-
-        Raises:
-            Warning: If interaction types are not found in the Hamiltonian parameters.
-        """
+    def _make_ham_inter(self):
         logger.debug(">>> Interaction._make_ham_inter")
 
         nx,ny,nz = self.lattice.shape
@@ -494,7 +395,7 @@ def _make_ham_inter(self):
         nd = norb * ns
 
         # Interaction Hamiltonian W[r,b,bp,a,ap]
-            #   H = W(r)^{\beta\beta^\prime\alpha\alpha^\prime}
+        #   H = W(r)^{\beta\beta^\prime\alpha\alpah^\prime}
         #        * c_{i\alpha}^\dagger c_{i\alpha^\prime} c_{j\beta^\prime}^\dagger c_{j\beta}
         ham_r = np.zeros((nx,ny,nz,*(ns,norb)*4), dtype=np.complex128)
 
@@ -506,17 +407,12 @@ def _make_ham_inter(self):
             'Ising':        { (0,0,0,0): 1, (1,1,1,1): 1, (0,0,1,1): -1, (1,1,0,0): -1 },
             'PairLift':     { (0,1,0,1): 1, (1,0,1,0): 1 },
             'Exchange':     { (0,1,1,0): -1, (1,0,0,1): -1 },
+            #--
             'PairHop':      { (0,0,1,1): 1, (1,1,0,0): 1 },
         }
 
         # coulomb-type interactions
         def _append_inter(type):
-            """
-            Append Coulomb-type interactions to the Hamiltonian.
-
-            Args:
-                type (str): Type of Coulomb interaction.
-            """
             logger.debug("_append_inter {}".format(type))
             spins = spin_table[type]
             for (irvec,orbvec), v in self.param_ham[type].items():
@@ -533,17 +429,11 @@ def _append_inter(type):
         #            c_{i\alpha\down}^\dagger c_{j\alpha^\prime\down}
         #        + (up <-> down)
         def _append_pairhop(type):
-            """
-            Append PairHop-type interactions to the Hamiltonian.
-
-            Args:
-                type (str): Type of PairHop interaction.
-            """
             spins = spin_table[type]
             for (irvec,orbvec), v in self.param_ham[type].items():
                 # take account of same-site interaction only
-                    if irvec == (0, 0, 0):
-                        a, b = orbvec
+                if (irvec == (0,0,0)):
+                    a, b = orbvec
                     for spinvec, w in spins.items():
                         s1,s2,s3,s4 = spinvec
                         # beta beta' alpha alpha'
@@ -606,14 +496,6 @@ class RPA:
     """
     @do_profile
     def __init__(self, param_ham, info_log, info_mode):
-        """
-        Initialize the RPA object.
-
-        Args:
-            param_ham (dict): Dictionary containing Hamiltonian parameters.
-            info_log (dict): Dictionary containing log information.
-            info_mode (dict): Dictionary containing mode information.
-        """
         logger.debug(">>> RPA.__init__")
 
         self.param_ham = param_ham
@@ -631,14 +513,7 @@ def __init__(self, param_ham, info_log, info_mode):
         pass
 
     def _set_scheme(self, info_mode):
-        """
-        Set the calculation scheme based on the provided mode information.
-
-        This method handles the `calc_scheme` parameter, which must be called after setting up interactions.
-
-        Args:
-            info_mode (dict): Dictionary containing mode information.
-        """
+        # handle calc_scheme: must be called after setting up interactions
 
         self.calc_scheme = info_mode.get("calc_scheme", "auto")
 
@@ -667,13 +542,6 @@ def _set_scheme(self, info_mode):
         self.enable_reduced = self.calc_scheme.lower() in ["reduced", "squashed"]
         
     def _init_param(self):
-        """
-        Initialize parameters for the RPA object.
-
-        This method sets up various parameters such as temperature, energy cutoff,
-        number of Matsubara frequencies, number of orbitals, spin degrees of freedom,
-        and coefficients for tail and external fields.
-        """
         logger.debug(">>> RPA._init_param")
 
         self.T = self.param_mod.get("T", 0.0)
@@ -750,16 +618,6 @@ def _init_param(self):
         pass
 
     def _round_to_int(self, val, mode):
-        """
-        Round a value to an integer based on the specified mode.
-
-        Args:
-            val (float): The value to be rounded.
-            mode (str): The rounding mode. Can be one of "as-is", "round", "round-up", or "round-down".
-
-        Returns:
-            int: The rounded integer value.
-        """
         import math
         mode = mode.lower()  # case-insensitive
         if   mode == "as-is":
@@ -791,9 +649,6 @@ def _round_to_int(self, val, mode):
         return ret
 
     def _show_params(self):
-        """
-        Log the RPA parameters.
-        """
         logger.debug(">>> RPA._show_params")
         logger.info("RPA parameters:")
         logger.info("    norbit          = {}".format(self.norb))
@@ -823,13 +678,6 @@ def _show_params(self):
 
     @do_profile
     def solve(self, green_info, path_to_output):
-        """
-        Solve the RPA calculations.
-
-        Args:
-            green_info (dict): Dictionary containing Green's function information.
-            path_to_output (str): Path to the output directory.
-        """
         logger.info("Start RPA calculations")
 
         beta = 1.0/self.T
@@ -1013,13 +861,6 @@ def solve(self, green_info, path_to_output):
 
     @do_profile
     def save_results(self, info_outputfile, green_info):
-        """
-        Save the RPA results to the specified output file.
-
-        Args:
-            info_outputfile (dict): Dictionary containing output file information.
-            green_info (dict): Dictionary containing Green's function information.
-        """
         logger.info("Save RPA results")
         path_to_output = info_outputfile["path_to_output"]
 
@@ -1051,15 +892,7 @@ def save_results(self, info_outputfile, green_info):
         pass
 
     def _init_wavevec(self):
-        """
-        Initialize wave vectors on sublatticed geometry.
-
-        This method calculates the reciprocal lattice vectors and the wave number table
-        for the sublatticed geometry.
-
-        The wave vectors are stored in `self.kvec` and the wave number table is stored
-        in `self.wavenum_table` and `self.wave_table`.
-        """
+        # wave vectors on sublatticed geometry
         def _klist(n):
             return np.roll( (np.arange(n)-(n//2)), -(n//2) )
 
@@ -1068,8 +901,7 @@ def _klist(n):
         omg = np.dot(rvec[0], np.cross(rvec[1], rvec[2]))
         kvec = np.array([
             np.cross(rvec[(i+1)%3], rvec[(i+2)%3])/omg * 2*np.pi/self.lattice.shape[i]
-            for i in range(3)
-        ])
+            for i in range(3) ])
 
         self.kvec = kvec  # store reciprocal lattice vectors
 
@@ -1085,23 +917,15 @@ def _klist(n):
                     v = kvec[0] * kx + kvec[1] * ky + kvec[2] * kz
                     wtable[ix,iy,iz] = v
         self.wave_table = wtable.reshape(nvol,3)
-    def _find_index_range(self, freq_range):
-        """
-        Decode Matsubara frequency index list.
-
-        This method decodes the Matsubara frequency index list based on the provided
-        frequency range. The frequency range can be specified as a single index, a
-        range with min and max, a range with min, max, and step, or a keyword.
-
-        Args:
-            freq_range (int, list, str): The frequency range specification.
 
-        Returns:
-            list: The list of frequency indices.
+    def _find_index_range(self, freq_range):
+        # decode matsubara frequency index list
+        #   1. single index
+        #   2. min, max, step
+        #   3. keyword
+        # note index n in [0 .. Nmat-1] corresponds to
+        #   w_n = (2*n-Nmat) * pi / beta
 
-        Raises:
-            ValueError: If the frequency range specification is invalid.
-        """
         nmat = self.nmat
 
         if type(freq_range) == int:
@@ -1134,22 +958,8 @@ def _find_index_range(self, freq_range):
 
         return freq_index
 
-    @do_profile
     @do_profile
     def read_init(self, info_inputfile):
-        """
-        Read initial configurations for the RPA calculations.
-
-        This method reads the initial configurations for the RPA calculations from
-        the specified input files. The configurations include chi0q, trans_mod, and
-        green_init.
-
-        Args:
-            info_inputfile (dict): Dictionary containing input file information.
-
-        Returns:
-            dict: Dictionary containing the read configurations.
-        """
         logger.info("RPA read initial configs")
         info = {}
 
@@ -1170,21 +980,6 @@ def read_init(self, info_inputfile):
         return info
 
     def _read_chi0q(self, file_name):
-        """
-        Read chi0q data from a file.
-
-        This method reads the chi0q data from the specified file and checks its size
-        and shape based on the calculation scheme.
-
-        Args:
-            file_name (str): Name of the file to read the chi0q data from.
-
-        Returns:
-            numpy.ndarray: The read chi0q data.
-
-        Raises:
-            SystemExit: If reading the file fails or the data shape is unexpected.
-        """
         logger.debug(">>> RPA._read_chi0q")
 
         try:
@@ -1265,15 +1060,6 @@ def _read_chi0q(self, file_name):
         return chi0q
 
     def _read_trans_mod(self, file_name):
-        """
-        Read the transformation module from a file and perform FFT.
-
-        Args:
-            file_name (str): Name of the file to read the transformation module from.
-
-        Returns:
-            numpy.ndarray: The transformed module in k-space.
-        """
         logger.debug(">>> RPA._read_trans_mod")
 
         try:
@@ -1298,15 +1084,6 @@ def _read_trans_mod(self, file_name):
         return tab_k
 
     def _read_green(self, file_name):
-        """
-        Read the Green's function from a file.
-
-        Args:
-            file_name (str): Name of the file to read the Green's function from.
-
-        Returns:
-            numpy.ndarray: The Green's function reshaped to the appropriate dimensions.
-        """
         logger.debug(">>> RPA._read_green")
         try:
             logger.info("read green from {}".format(file_name))
@@ -1320,17 +1097,10 @@ def _read_green(self, file_name):
         nvol = self.lattice.nvol
         nd = self.nd
         return green.reshape(nvol,nd,nd)
-    def _reshape_green(self, green_):
-        def _reshape_green(self, green_):
-            """
-            Convert the Green's function into sublattice format.
 
-            Args:
-                green_ (numpy.ndarray): The original Green's function.
+    def _reshape_green(self, green_):
+        # convert green function into sublattice
 
-            Returns:
-                numpy.ndarray: The Green's function reshaped for the sublattice.
-            """
         Lx,Ly,Lz = self.lattice.cellshape
         Lvol = Lx * Ly * Lz
         Bx,By,Bz = self.lattice.subshape
@@ -1383,15 +1153,6 @@ def _unpack_index(x, n):
         return green_sub
 
     def _calc_trans_mod(self, g0):
-        """
-        Calculate the transformation module.
-
-        Args:
-            g0 (numpy.ndarray): The Green's function.
-
-        Returns:
-            numpy.ndarray: The transformed Hamiltonian in k-space.
-        """
         logger.debug(">>> RPA._calc_trans_mod")
 
         nx,ny,nz = self.lattice.shape
@@ -1419,19 +1180,6 @@ def _calc_trans_mod(self, g0):
 
     @do_profile
     def _calc_epsilon_k(self, green_info):
-        """
-        Calculate the epsilon\_k values for the RPA calculations.
-
-        This method calculates the epsilon\_k values based on the provided Green's function information.
-        It determines the transfer term H0(k) and performs necessary transformations based on the spin-orbital
-        interactions and external fields.
-
-        Args:
-            green_info (dict): Dictionary containing Green's function information.
-
-        Returns:
-            None
-        """
         logger.debug(">>> RPA._calc_epsilon_k")
 
         nvol = self.lattice.nvol
@@ -1514,20 +1262,6 @@ def _calc_epsilon_k(self, green_info):
 
     @do_profile
     def _find_mu(self, Ncond, T):
-        """
-        Find the chemical potential (mu) for the given number of electrons and temperature.
-
-        This method calculates the chemical potential (mu) such that the number of electrons
-        matches the specified value (Ncond) at the given temperature (T). It uses the eigenvalues
-        and eigenvectors of the Hamiltonian to perform the calculation.
-
-        Args:
-            Ncond (float): The number of electrons.
-            T (float): The temperature.
-
-        Returns:
-            tuple: A tuple containing the distribution (dist) and the chemical potential (mu).
-        """
         logger.debug(">>> RPA._find_mu")
 
         from scipy import optimize
@@ -1544,17 +1278,6 @@ def _find_mu(self, Ncond, T):
         occupied_number = Ncond
 
         def _fermi(t, mu, ev):
-            """
-            Calculate the Fermi-Dirac distribution.
-
-            Args:
-                t (float): Temperature.
-                mu (float): Chemical potential.
-                ev (numpy.ndarray): Eigenvalues.
-
-            Returns:
-                numpy.ndarray: Fermi-Dirac distribution values.
-            """
             w_ = (ev - mu) / t
             mask_ = w_ < ene_cutoff
             w1_ = np.where( mask_, w_, 0.0 )
@@ -1563,20 +1286,11 @@ def _fermi(t, mu, ev):
             return v_
 
         def _calc_delta_n(mu):
-            """
-            Calculate the difference between the number of electrons and the target number.
-
-            Args:
-                mu (float): Chemical potential.
-
-            Returns:
-                float: Difference between the calculated and target number of electrons.
-            """
             ff = _fermi(T, mu, w)
             nn = np.einsum('gkal,gkl,gkal->', np.conjugate(v), ff, v)
             return nn.real - occupied_number
 
-        # Find mu such that <n>(mu) = N0
+        # find mu s.t. <n>(mu) = N0
         is_converged = False
         if (_calc_delta_n(ev[0]) * _calc_delta_n(ev[-1])) < 0.0:
             logger.debug("RPA._find_mu: try bisection")
@@ -1599,14 +1313,10 @@ def _calc_delta_n(mu):
     @do_profile
     def _calc_green(self, beta, mu):
         """
-        Calculate the Green's function for the RPA calculations.
-
-        Args:
-            beta (float): Inverse temperature.
-            mu (float): Chemical potential.
-
-        Returns:
-            tuple: A tuple containing the Green's function and the tail correction.
+        ew: eigenvalues  ew[g,k,i]    i-th eigenvalue of wavenumber k in block g
+        ev: eigenvectors ev[g,k,a,i]  i-th eigenvector corresponding to ew[g,k,i]
+        beta: inverse temperature
+        mu: chemical potential
         """
         logger.debug(">>> RPA._calc_green")
 
@@ -1637,19 +1347,15 @@ def _calc_green(self, beta, mu):
         green_tail = np.einsum('gkaj,gkbj,gl->glkab', ev, np.conj(ev), np.tile(np.ones(nmat), (nblock,1))) * aa * 0.5 * beta
 
         return green, green_tail
+
     @do_profile
     def _calc_chi0q(self, green_kw, green0_tail, beta):
-        """
-        Calculate the chi0q values for the RPA calculations.
-
-        Args:
-            green_kw (numpy.ndarray): Green's function in Matsubara frequency space.
-            green0_tail (numpy.ndarray): Tail correction for the Green's function.
-            beta (float): Inverse temperature.
+        # green function
+        #   green_kw[g,l,k,a,b]: G_ab(k,ie) in block g
+        #     l : matsubara freq i\epsilon_l = i\pi(2*l+1-N)/\beta for l=0..N-1
+        #     k : wave number kx,ky,kz
+        #     a,b : orbital and spin index a=(s,alpha) b=(t,beta)
 
-        Returns:
-            numpy.ndarray: The chi0q values in Matsubara frequency space.
-        """
         logger.debug(">>> RPA._calc_chi0q")
 
         nx,ny,nz = self.lattice.shape
@@ -1703,18 +1409,9 @@ def _calc_chi0q(self, green_kw, green0_tail, beta):
             axis=1).reshape(nblock,nmat,nvol,*nd_shape) * (-1.0/beta)
 
         return chi0_qw
+
     @do_profile
     def _solve_rpa(self, chi0q, ham):
-        """
-        Solve the RPA equations.
-
-        Args:
-            chi0q (numpy.ndarray): The chi0q values.
-            ham (numpy.ndarray): The Hamiltonian in k-space.
-
-        Returns:
-            numpy.ndarray: The solution to the RPA equations.
-        """
         logger.debug(">>> RPA._solve_rpa")
 
         nvol = self.lattice.nvol
@@ -1738,16 +1435,6 @@ def _solve_rpa(self, chi0q, ham):
 
 
 def run(*, input_dict: Optional[dict] = None, input_file: Optional[str] = None):
-    """
-    Run the RPA calculations based on the provided input dictionary and file.
-
-    Args:
-        input_dict (Optional[dict]): Dictionary containing input parameters.
-        input_file (Optional[str]): Path to the input file (not used in the current implementation).
-
-    Raises:
-        RuntimeError: If input_dict is not provided.
-    """
     logger = logging.getLogger("run")
 
     if input_dict is None: