class Approach(object):
    """
    Sample class used to define different approaches.

    Attributes
    ----------
    qd : QuantumDot
        QuantumDot object.
    leads : LeadsTunneling
        LeadsTunneling object.
    si : StateIndexingDM
        StateIndexingDM object.
    kern : array
        Kernel (Liouvillian) representing the master equation.
    bvec : array
        Right hand side column vector for master equation.
        The entry funcp.norm_row is 1 representing normalization condition.
    sol0 : array
        Least squares solution for the master equation.
    phi0 : array
        Values of zeroth order density matrix elements.
    phi1 : array
        Values of first order density matrix elements
        stored in nleads by ndm1 numpy array.
    current : array
        Values of the current having nleads entries.
    energy_current : array
        Values of the energy current having nleads entries.
    heat_current : array
        Values of the heat current having nleads entries.
    funcp : FunctionProperties
        FunctionProperties object.
    paulifct : array
        Factors used for generating Pauli master equation kernel.
    phi1fct : array
        Factors used for generating 1vN, Redfield master equation kernels.
    phi1fct_energy : array
        Factors used to calculate energy and heat currents in 1vN, Redfield approaches.
    tLba : array
        Jump operator matrix in many-body basis for Lindblad approach.
    """

    kerntype = 'not defined'
    dtype = doublenp
    indexing_class_name = 'StateIndexingDM'

    #region Properties

    @property
    def solmethod(self):
        return self.funcp.solmethod
    @solmethod.setter
    def solmethod(self, value):
        self.funcp.solmethod = value
        self.prepare_solver()

    @property
    def mfreeq(self):
        return self.funcp.mfreeq
    @mfreeq.setter
    def mfreeq(self, value):
        self.funcp.mfreeq = value

    @property
    def symq(self):
        return self.funcp.symq
    @symq.setter
    def symq(self, value):
        if value == self.funcp.symq:
            return
        self.is_prepared = False
        self.funcp.symq = value

    @property
    def norm_row(self):
        return self.funcp.norm_row
    @norm_row.setter
    def norm_row(self, value):
        self.funcp.norm_row = value

    @property
    def itype(self):
        return self.funcp.itype
    @itype.setter
    def itype(self, value):
        self.funcp.itype = value

    #endregion Properties

    def __init__(self, builder):
        """
        Initialization of the Approach class.

        Parameters
        ----------
        builder : Builder, Approach, etc. object.
            Any object having qd, leads, si, and funcp attributes.
        """
        self.qd = builder.qd
        self.leads = builder.leads
        self.si = builder.si
        self.funcp = builder.funcp
        self.restart()

    def restart(self):
        """Restart values of some variables."""
        self.is_prepared = False
        self.make_kern_copy = False

        self.kern, self.bvec, self.norm_vec = None, None, None
        self.sol0, self.phi0, self.phi1 = None, None, None
        self.dphi0_dt = None
        self.current = None
        self.energy_current = None
        self.heat_current = None
        self.phi1fct, self.paulifct = None, None
        self.phi1fct_energy = None
        self.tLba = None
        self.kernel_handler = None

    #region Preparation

    def get_kern_size(self):
        return self.si.ndm0r

    def prepare_kern(self):
        if self.is_prepared and not self.si.states_changed:
            self.clean_arrays()
            return

        self.prepare_kernel_handler()
        self.prepare_arrays()
        self.prepare_solver()

        self.si.states_changed = False
        self.is_prepared = True

    def clean_arrays(self):
        if not self.mfreeq:
            self.kern.fill(0.0)
            self.bvec.fill(0.0)

        self.current.fill(0.0)
        self.energy_current.fill(0.0)
        self.heat_current.fill(0.0)

    def prepare_arrays(self):
        kern_size = self.get_kern_size()
        nleads = self.si.nleads

        self.current = np.zeros(nleads, dtype=doublenp)
        self.energy_current = np.zeros(nleads, dtype=doublenp)
        self.heat_current = np.zeros(nleads, dtype=doublenp)

        self.phi0 = np.zeros(kern_size, dtype=self.dtype)
        self.kernel_handler.set_phi0(self.phi0)

        if self.funcp.mfreeq:
            self.dphi0_dt = np.zeros(kern_size, dtype=self.dtype)
            self.kernel_handler.set_dphi0_dt(self.dphi0_dt)
        else:
            kern_size_rows = kern_size if self.funcp.symq else kern_size+1

            self.kern = np.zeros((kern_size_rows, kern_size), dtype=self.dtype, order='F')

            self.kernel_handler.set_kern(self.kern)

            self.bvec = np.zeros(kern_size_rows, dtype=self.dtype)

            self.replaced_eq = np.zeros(kern_size, dtype=self.dtype)

            self.norm_vec = np.zeros(kern_size, dtype=self.dtype)

            self.generate_norm_vec()

    def prepare_kernel_handler(self):
        if self.funcp.mfreeq:
            self.kernel_handler = KernelHandlerMatrixFree(self.si)
        else:
            self.kernel_handler = KernelHandler(self.si)

    def prepare_solver(self):
        solmethod = self.funcp.solmethod
        # Determine the proper solution method
        if self.funcp.mfreeq:
            solmethod = solmethod if solmethod is not None else 'krylov'
            self.funcp.solmethod = solmethod
            return

        symq = self.funcp.symq
        if solmethod is None:
            solmethod = 'solve' if symq else 'lsqr'
        if not symq and solmethod != 'lsqr' and solmethod != 'lsmr':
            print("WARNING: Using solmethod=lsqr, because the kernel is not symmetric, symq=False.")
            solmethod = 'lsqr'

        self.funcp.solmethod = solmethod

    #endregion Preparation

    #region Generation

    def generate_norm_vec(self):
        """
        Generates normalisation condition for 1vN approach.

        Parameters
        ----------
        norm_vec : array
            (Modifies) Left hand side of the normalisation condition.
        """
        si = self.si

        for charge in range(si.ncharge):
            for b in si.statesdm[charge]:
                bb = si.get_ind_dm0(b, b, charge)
                self.norm_vec[bb] += 1

    def generate_fct(self):
        pass

    def generate_kern(self):
        """
        Generates a kernel (Liouvillian) matrix corresponding to first order von Neumann approach (1vN).

        Parameters
        ----------
        kern : array
            (Modifies) Kernel matrix for 1vN approach.
        """
        E = self.qd.Ea
        si, kh = self.si, self.kernel_handler
        ncharge, statesdm = si.ncharge, si.statesdm

        for bcharge in range(ncharge):
            for b, bp in itertools.combinations_with_replacement(statesdm[bcharge], 2):
                if not (kh.is_included(b, bp, bcharge) and kh.is_unique(b, bp, bcharge)):
                    continue
                kh.set_energy(E[b]-E[bp], b, bp, bcharge)
                self.generate_coupling_terms(b, bp, bcharge)

    def generate_coupling_terms(self, b, bp, bcharge):
        pass

    def generate_current(self):
        pass

    def generate_vec(self, phi0):
        """
        Acts on given phi0 with Liouvillian of first-order approach.

        Parameters
        ----------
        phi0 : array
            Some values of zeroth order density matrix elements.
        self : Approach
            Approach object.

        Returns
        -------
        dphi0_dt : array
            Values of zeroth order density matrix elements
            after acting with Liouvillian, i.e., dphi0_dt=L(phi0p).
        """
        norm_row = self.funcp.norm_row

        kh = self.kernel_handler
        kh.set_phi0(phi0)
        norm = kh.get_phi0_norm()

        self.dphi0_dt.fill(0.0)

        # Here dphi0_dt and norm will be implicitly calculated by using KernelHandlerMatrixFree
        self.generate_kern()

        self.dphi0_dt[norm_row] = norm-1

        return self.dphi0_dt

    #endregion Generation

    #region Solution

    def solve_kern(self):
        """Finds the stationary state using least squares or using LU decomposition."""
        solmethod = self.funcp.solmethod
        symq = self.funcp.symq
        norm_row = self.funcp.norm_row
        replaced_eq = self.replaced_eq

        kern = self.kern
        bvec = self.bvec

        # Replace one equation by the normalisation condition
        if symq:
            replaced_eq[:] = kern[norm_row]
            kern[norm_row] = self.norm_vec
            bvec[norm_row] = 1
        else:
            kern[-1] = self.norm_vec
            bvec[-1] = 1

        # Try to solve the master equation
        try:
            if solmethod == 'solve':
                self.sol0 = [np.linalg.solve(kern, bvec)]
            elif solmethod == 'lsqr':
                self.sol0 = np.linalg.lstsq(kern, bvec, rcond=-1)

            self.phi0[:] = self.sol0[0]
            self.success = True
        except Exception as exept:
            self.funcp.print_error(exept)
            self.phi0.fill(0.0)
            self.success = False

        # Return back the replaced equation
        if symq:
            kern[norm_row] = replaced_eq

    def solve_matrix_free(self):
        """Finds the stationary state using matrix free methods like broyden, krylov, etc."""
        solmethod = self.funcp.solmethod
        phi0_init = self.funcp.phi0_init

        if phi0_init is None:
            self.funcp.print_warning(0, "WARNING: For mfreeq=True no phi0_init is specified. " +
                                     "Using phi0_init[0]=1.0 as a default. " +
                                     "This warning will not be shown again.")
            phi0_init = self.set_phi0_init()

        try:
            self.sol0 = optimize.root(self.generate_vec, phi0_init, method=solmethod)
            self.phi0[:] = self.sol0.x
            self.success = self.sol0.success
        except Exception as exept:
            self.funcp.print_error(exept)
            self.phi0.fill(0.0)
            self.success = False

    def set_phi0_init(self):
        phi0_init = np.zeros(self.get_kern_size(), dtype=self.dtype)
        phi0_init[0] = 1.0
        return phi0_init

    def rotate(self):
        self.leads.rotate(self.qd.vecslst)

    def solve(self, qdq=True, rotateq=True, masterq=True, currentq=True, *args, **kwargs):
        """
        Solves the master equation.

        Parameters
        ----------
        qdq : bool
            Diagonalise many-body quantum dot Hamiltonian
            and express the lead matrix Tba in the eigenbasis.
        rotateq : bool
            Rotate the many-body tunneling matrix Tba.
        masterq : bool
            Solve the master equation.
        currentq : bool
            Calculate the current.
        """
        if qdq:
            self.qd.diagonalise()
            if rotateq:
                self.rotate()
        #
        if masterq:
            self.prepare_kern()
            self.generate_fct()
            if not self.funcp.mfreeq:
                self.generate_kern()
                self.solve_kern()
            else:
                self.solve_matrix_free()
            if currentq:
                self.generate_current()

    #endregion Solution class StateIndexingDM(StateIndexing):
    """
    Class for indexing density matrix elements.
    Separates the off-diagonal density matrix elements into real and imaginary parts.
    Derived from StateIndexing.

    Attributes
    ----------
    ndm0 : int
        Number of density matrix elements to zeroth order corresponding to statesdm,
        when symmetries are used.
    ndm0r : int
        Number of unique density matrix elements to zeroth order corresponding to statesdm,
        when symmetries are used.
        The real and imaginary parts of the off-diagonal matrix element are considered
        as separate entities.
    ndm0_ : int
        Number of density matrix elements to zeroth order corresponding to statesdm,
        when no symmetries are used.
    ndm0_tot : int
        Total number of density matrix elements to zeroth order for given nsingle,
        when no symmetries are used.
    ndm1, ndm1_ : int, int
        Number of density matrix elements to first order corresponding to statesdm,
        when no symmetries are used.
    ndm1_tot : int
        Total number of density matrix elements to first order for given nsingle,
        when no symmetries are used.
    npauli_ : int
        Number of diagonal density matrix elements to zeroth order corresponding to statesdm,
        when no symmetries are used.
    npauli : int
        Number of diagonal density matrix elements to zeroth order corresponding to statesdm,
        when symmetries are used.
    statesdm : list
        List of states considered in calculations involving density matrix.
    shiftlst0 : numpy array
        Array necessary for labeling zeroth order density matrix elements. See method get_ind_dm0().
    shiftlst1 : numpy array
        Array necessary for labeling first order density matrix elements. See method get_ind_dm1().
    lenlst : numpy array
        Array containing number of many-body states for each charge. lenlst[charge]=number of states.
    dictdm : numpy array
        Array giving enumeration of considered many-body states in statesdm.
    inddm0 : dict
        Dictionary between index of density matrix element and pair of many-body states (a,b).
    mapdm0 : list
        List showing which density matrix elements are mapped to each other due to symmetries.
        For example mapdm0[1]=1, mapdm0[2]=1, shows that 1 and 2 are equivalent.
        Also mapdm0[ind]=-1 shows that the state is excluded.
    booldm0 : list
        List giving density matrix elements which will be used.
        Other elements will be mapped to elements, which have booldm0[ind]=True.
        From example for mapdm0 we have booldm0[1]=True, booldm0[2]=False.
    conjdm0 : list
        List showing, which density matrix elements are complex conjugate and are not unique.
    """

    def __init__(self, nsingle, indexing='Lin', symmetry=None, nleads=0):
        """
        Initialization of the StateIndexingDM class

        Parameters
        ----------
        nsingle : int
            Number of single-particle states.
        indexing : str
            String determining type of the indexing. Possible values are 'Lin', 'charge', 'sz', 'ssq'.
        symmetry : str
            String determining that the states will be augmented by the symmetry.
        """
        StateIndexing.__init__(self, nsingle, indexing, symmetry, nleads)
        self.ndm0_tot = int(factorial(2*self.nsingle)/factorial(self.nsingle)**2)
        self.ndm1_tot = int(self.nsingle/(self.nsingle+1)*self.ndm0_tot)
        #
        self.shiftlst0 = np.zeros(self.ncharge+1, dtype=longnp)
        self.shiftlst1 = np.zeros(self.ncharge, dtype=longnp)
        self.lenlst = np.zeros(self.ncharge, dtype=longnp)
        self.dictdm = np.zeros(self.nmany, dtype=longnp)
        #
        self.statesdm, self.mapdm0, self.mapdm1 = None, None, None
        self.set_statesdm(self.chargelst)

    def set_statesdm(self, statesdm):
        """
        Sets the considered many-body states for density matrix calculations.

        Parameters
        ----------
        statesdm: list
            List containing indices of many-body state under consideration.
        """
        self.statesdm = list(statesdm)
        self.statesdm.append([])
        self.ndm0_, self.ndm1_, self.npauli_ = 0, 0, 0
        for j1 in range(self.ncharge):
            self.npauli_ += len(statesdm[j1])
            self.ndm0_ += len(statesdm[j1])**2
            if j1 < self.ncharge-1:
                self.ndm1_ += len(statesdm[j1])*len(statesdm[j1+1])
        self.set_dictdm()

        self.states_changed = True

    def set_dictdm(self):
        """
        Makes dictdm, shiftlst0, and shiftlst1 necessary for density-matrix element indexing.
        """
        for j1 in range(self.ncharge):
            self.lenlst[j1] = len(self.statesdm[j1])
            self.shiftlst0[j1+1] = self.shiftlst0[j1] + len(self.statesdm[j1])**2
            if j1 < self.ncharge-1:
                self.shiftlst1[j1+1] = (self.shiftlst1[j1]
                                        + len(self.statesdm[j1])*len(self.statesdm[j1+1]))
            counter = 0
            for j2 in self.statesdm[j1]:
                self.dictdm[j2] = counter
                counter += 1
        self.set_mapdm()

    def set_mapdm(self):
        """
        Reduce the number of diagonal matrix elements by using symmetries.
        """
        # noinspection PyShadowingNames
        def add_elem(counter, b, bp, charge, dictq=True, conjq=True):
            bbp = self.get_ind_dm0(b, bp, charge, maptype=0)
            self.mapdm0[bbp] = counter
            self.booldm0[bbp] = dictq
            self.conjdm0[bbp] = conjq
            if dictq:
                self.inddm0.update({counter: (b, bp)})
            return counter+1
        #
        self.mapdm0 = np.ones(self.ndm0_, dtype=longnp)*(-1)
        self.booldm0 = np.zeros(self.ndm0_, dtype=boolnp)
        self.conjdm0 = np.zeros(self.ndm0_, dtype=boolnp)
        self.inddm0 = {}
        counter = 0
        # Diagonal density matrix elements
        for charge in range(self.ncharge):
            for b in self.statesdm[charge]:
                if self.indexing == 'ssq':
                    (b_ch,  b_sz,  b_ssq,  b_alpha) = self.ind_qn[b]
                    if b_ssq == -b_sz:
                        add_elem(counter, b, b, charge)
                        for sz in range(b_sz+2, b_ssq+1, 2):
                            b1 = self.qn_ind[(b_ch,  sz, b_ssq,  b_alpha)]
                            add_elem(counter, b1, b1, charge, dictq=False)
                        counter = counter+1
                else:
                    counter = add_elem(counter, b, b, charge)
        self.npauli = counter
        # Off-diagonal density matrix elements
        for charge in range(self.ncharge):
            for b, bp in itertools.combinations(self.statesdm[charge], 2):
                if self.indexing == 'sz':
                    if self.ind_qn[b][1] == self.ind_qn[bp][1]:
                        add_elem(counter, bp, b, charge, dictq=False, conjq=False)
                        counter = add_elem(counter, b, bp, charge)
                elif self.indexing == 'ssq':
                    (b_ch,  b_sz,  b_ssq,  b_alpha) = self.ind_qn[b]
                    (bp_ch, bp_sz, bp_ssq, bp_alpha) = self.ind_qn[bp]
                    if (b_sz, b_ssq) == (bp_sz, bp_ssq):
                        if b_ssq == -b_sz:
                            add_elem(counter, bp, b, charge, dictq=False, conjq=False)
                            add_elem(counter, b, bp, charge)
                            for sz in range(b_sz+2, b_ssq+1, 2):
                                b1 = self.qn_ind[(b_ch,  sz, b_ssq,  b_alpha)]
                                b1p = self.qn_ind[(bp_ch, sz, bp_ssq, bp_alpha)]
                                add_elem(counter, b1p, b1, charge, dictq=False, conjq=False)
                                add_elem(counter, b1, b1p, charge, dictq=False)
                            counter = counter+1
                else:
                    add_elem(counter, bp, b, charge, dictq=False, conjq=False)
                    counter = add_elem(counter, b, bp, charge)
        self.ndm0 = counter
        self.ndm0r = self.npauli+2*(self.ndm0-self.npauli)
        self.ndm1 = self.ndm1_

    def get_ind_dm0(self, b, bp, charge, maptype=1):
        """
        Get the index of zeroth order density matrix element.

        Parameters
        ----------
        b,bp : int
            Indices of many-body states within the same charge state.
        charge : int
            Charge of b and bp.
        maptype : int
            Determines what kind of mapping to use when returning the index.

        Returns
        -------
        int
            Index of the zeroth order density matrix element.
        """
        # l = len(self.statesdm[charge])
        # i = self.dictdm[b]
        # j = self.dictdm[bp]
        # shiftas = self.shiftlst0[charge]
        # print('index is', l*i + j + shiftas)
        if maptype == 0:
            return (self.lenlst[charge]*self.dictdm[b] + self.dictdm[bp]
                    + self.shiftlst0[charge])
        elif maptype == 1:
            return (self.mapdm0[self.lenlst[charge]*self.dictdm[b] + self.dictdm[bp]
                    + self.shiftlst0[charge]])
        elif maptype == 2:
            return (self.booldm0[self.lenlst[charge]*self.dictdm[b] + self.dictdm[bp]
                    + self.shiftlst0[charge]])
        elif maptype == 3:
            return (self.conjdm0[self.lenlst[charge]*self.dictdm[b] + self.dictdm[bp]
                    + self.shiftlst0[charge]])

    def get_ind_dm1(self, c, b, bcharge):
        """
        Get the index of first order density matrix element.

        Parameters
        ----------
        c,b : int
            Indices of many-body states differing by single charge, charge(c)>charge(b).
        bcharge : int
            Charge of many-body state b.

        Returns
        -------
        int
            Index of the first order density matrix element.
        """
        # l = len(self.statesdm[charge])
        # i = self.dictdm[c]
        # j = self.dictdm[b]
        # shiftas = self.shiftlst1[charge]
        # print('index is', l*i + j + shiftas)
        #       b
        #     -----
        #   c |   |
        #     -----
        return self.lenlst[bcharge]*self.dictdm[c] + self.dictdm[b] + self.shiftlst1[bcharge]

    def remove_fock_states(self, lin_state_indices):
        StateIndexing.remove_fock_states(self, lin_state_indices)
        self.set_statesdm(self.chargelst) class StateIndexingDMc(StateIndexing):
    """
    Class for indexing density matrix elements.
    Does not separates the off-diagonal density matrix elements into real and imaginary parts,
    and treats Phi[0]_{b,bp} and Phi[0]_{bp,b} as separate entities.
    Derived from StateIndexing.

    Attributes
    ----------
    Same as in StateIndexingDM.
    """

    def __init__(self, nsingle, indexing='Lin', symmetry=None, nleads=0):
        """
        Initialization of the StateIndexingDMc class

        Parameters
        ----------
        nsingle : int
            Number of single-particle states.
        indexing : str
            String determining type of the indexing. Possible values are 'Lin', 'charge', 'sz', 'ssq'.
        symmetry : str
            String determining that the states will be augmented by the symmetry.
        """
        StateIndexing.__init__(self, nsingle, indexing, symmetry, nleads)
        self.ndm0_tot = int(factorial(2*self.nsingle)/factorial(self.nsingle)**2)
        self.ndm1_tot = int(self.nsingle/(self.nsingle+1)*self.ndm0_tot)
        #
        self.shiftlst0 = np.zeros(self.ncharge+1, dtype=longnp)
        self.shiftlst1 = np.zeros(self.ncharge, dtype=longnp)
        self.lenlst = np.zeros(self.ncharge, dtype=longnp)
        self.dictdm = np.zeros(self.nmany, dtype=longnp)
        #
        self.statesdm, self.mapdm0, self.mapdm1 = None, None, None
        self.set_statesdm(self.chargelst)

    def set_statesdm(self, statesdm):
        """
        Sets the considered many-body states for density matrix calculations.

        Parameters
        ----------
        statesdm: list
            List containing indices of many-body state under consideration.
        """
        self.statesdm = statesdm
        self.statesdm.append([])
        self.ndm0_, self.ndm1_, self.npauli_ = 0, 0, 0
        for j1 in range(self.ncharge):
            self.npauli_ += len(statesdm[j1])
            self.ndm0_ += len(statesdm[j1])**2
            if j1 < self.ncharge-1:
                self.ndm1_ += len(statesdm[j1])*len(statesdm[j1+1])
        self.set_dictdm()

        self.states_changed = True

    def set_dictdm(self):
        """
        Makes dictdm, shiftlst0, and shiftlst1 necessary for density-matrix element indexing.
        """
        for j1 in range(self.ncharge):
            self.lenlst[j1] = len(self.statesdm[j1])
            self.shiftlst0[j1+1] = self.shiftlst0[j1] + len(self.statesdm[j1])**2
            if j1 < self.ncharge-1:
                self.shiftlst1[j1+1] = (self.shiftlst1[j1]
                                        + len(self.statesdm[j1])*len(self.statesdm[j1+1]))
            counter = 0
            for j2 in self.statesdm[j1]:
                self.dictdm[j2] = counter
                counter += 1
        self.set_mapdm()

    def set_mapdm(self):
        """
        Reduce the number of diagonal matrix elements by using symmetries.
        """
        # noinspection PyShadowingNames
        def add_elem(counter, b, bp, charge, dictq=True):
            bbp = self.get_ind_dm0(b, bp, charge, maptype=0)
            self.mapdm0[bbp] = counter
            self.booldm0[bbp] = dictq
            if dictq:
                self.inddm0.update({counter: (b, bp)})
            return counter+1
        #
        self.mapdm0 = np.ones(self.ndm0_, dtype=longnp)*(-1)
        self.booldm0 = np.zeros(self.ndm0_, dtype=boolnp)
        self.conjdm0 = None
        self.inddm0 = {}
        counter = 0
        # Diagonal density matrix elements
        for charge in range(self.ncharge):
            for b in self.statesdm[charge]:
                if self.indexing == 'ssq':
                    (b_ch,  b_sz,  b_ssq,  b_alpha) = self.ind_qn[b]
                    if b_ssq == -b_sz:
                        add_elem(counter, b, b, charge)
                        for sz in range(b_sz+2, b_ssq+1, 2):
                            b1 = self.qn_ind[(b_ch,  sz, b_ssq,  b_alpha)]
                            add_elem(counter, b1, b1, charge, dictq=False)
                        counter = counter+1
                else:
                    counter = add_elem(counter, b, b, charge)
        self.npauli = counter
        # self.npauli = 0
        # Off-diagonal density matrix elements
        for charge in range(self.ncharge):
            # for b, bp in itertools.product(self.statesdm[charge], self.statesdm[charge]):
            for b, bp in itertools.permutations(self.statesdm[charge], 2):
                if self.indexing == 'sz':
                    if self.ind_qn[b][1] == self.ind_qn[bp][1]:
                        counter = add_elem(counter, b, bp, charge)
                elif self.indexing == 'ssq':
                    (b_ch,  b_sz,  b_ssq,  b_alpha) = self.ind_qn[b]
                    (bp_ch, bp_sz, bp_ssq, bp_alpha) = self.ind_qn[bp]
                    if (b_sz, b_ssq) == (bp_sz, bp_ssq):
                        if b_ssq == -b_sz:
                            add_elem(counter, b, bp, charge)
                            for sz in range(b_sz+2, b_ssq+1, 2):
                                b1 = self.qn_ind[(b_ch,  sz, b_ssq,  b_alpha)]
                                b1p = self.qn_ind[(bp_ch, sz, bp_ssq, bp_alpha)]
                                add_elem(counter, b1, b1p, charge, dictq=False)
                            counter = counter+1
                else:
                    counter = add_elem(counter, b, bp, charge)
        self.ndm0 = counter
        self.ndm0r = self.npauli+2*(self.ndm0-self.npauli)
        self.ndm1 = self.ndm1_

    def get_ind_dm0(self, b, bp, charge, maptype=1):
        """
        Get the index of zeroth order density matrix element.

        Parameters
        ----------
        b,bp : int
            Indices of many-body states within the same charge state.
        charge : int
            Charge of b and bp.
        maptype : int
            Determines what kind of mapping to use when returning the index.

        Returns
        -------
        int
            Index of the zeroth order density matrix element.
        """
        # l = len(self.statesdm[charge])
        # i = self.dictdm[b]
        # j = self.dictdm[bp]
        # shiftas = self.shiftlst0[charge]
        # print('index is', l*i + j + shiftas)
        if maptype == 0:
            return self.lenlst[charge]*self.dictdm[b] + self.dictdm[bp] + self.shiftlst0[charge]
        elif maptype == 1:
            return (self.mapdm0[self.lenlst[charge]*self.dictdm[b] + self.dictdm[bp]
                    + self.shiftlst0[charge]])
        elif maptype == 2:
            return (self.booldm0[self.lenlst[charge]*self.dictdm[b] + self.dictdm[bp]
                    + self.shiftlst0[charge]])

    def get_ind_dm1(self, c, b, bcharge):
        """
        Get the index of first order density matrix element.

        Parameters
        ----------
        c,b : int
            Indices of many-body states differing by single charge, charge(c)>charge(b).
        bcharge : int
            Charge of many-body state b.

        Returns
        -------
        int
            Index of the first order density matrix element.
        """
        # l = len(self.statesdm[charge])
        # i = self.dictdm[c]
        # j = self.dictdm[b]
        # shiftas = self.shiftlst1[charge]
        # print('index is', l*i + j + shiftas)
        #       b
        #     -----
        #   c |   |
        #     -----
        return self.lenlst[bcharge]*self.dictdm[c] + self.dictdm[b] + self.shiftlst1[bcharge]

    def remove_fock_states(self, lin_state_indices):
        StateIndexing.remove_fock_states(self, lin_state_indices)
        self.set_statesdm(self.chargelst) class QuantumDot(object):
    """
    Class for constructing and diagonalising many-body Hamiltonian describing
    the quantum dot.

    Attributes
    ----------
    hsingle : dict, list, or ndarray
        Dictionary, list, or array corresponding to single-particle hopping (tunneling) Hamiltonian.
        On input list or array gets converted to dictionary.
    coulomb : dict, list, or ndarray
        Dictionary, list, or array or containing coulomb matrix elements.
        For dictionary:    coulomb[(m, n, k, l)] = U, where m, n, k, l are the state labels.
        For list or array: coulomb[i] is list of the format [m, n, k, l, U].
        U is the strength of the coulomb interaction between the states (m, n, k, l).
        Note that only the matrix elements k>l, n>m have to be specified.
        On input list or array gets converted to dictionary.
    si : StateIndexing
        StateIndexing or StateIndexingDM object.
    herm_hs : bool
        When herm_hs=True the conjugated elements are added automatically to many-body Hamiltonian
        from hsingle. For example, when hsingle={(0,1): val} the element {(1,0): val.conjugate()}
        will also be added.
    herm_c : bool
        When herm_c=True the conjugated elements are added automatically to many-body Hamiltonian
        from coulomb. For example, when coulomb={(0,2,3,1): val} the element
        {(1,3,2,0): val.conjugate()} will also be added. However, we note that for m_less_n=True
        if coulomb={(0,2,1,3): val} the element (3,1,2,0) will not be included.
    m_less_n : bool
        When m_less_n=True the coulomb matrix element (m, n, k, l) has to have m<n.
    mtype : type
        Defines the type of the Hamiltonian matrices. For example, float, complex, etc.
    ham_ssq_q : bool
        If ham_ssq_q is False then for 'ssq' indexing only Hamiltonian with lowest abs(sz) value
        is stored.
    valslst : list of ndarrays
        List containing ncharge entries, where valslst[charge] is an array having Hamiltonian
        eigenvalues for definite charge.
    vecslst : list of ndarrays
        List containing ncharge entries, where vecslst[charge] is an array corresponding to
        Hamiltonian eigenvectors for definite charge.
    hamlst : list of ndarrays
        List containing ncharge entries, where hamlst[charge] is an array corresponding to many-body
        Hamiltonian for definite charge.
    Ea : ndarray
        nmany by 1 array containing many-body Hamiltonian eigenvalues.
    """

    def __init__(self, hsingle, coulomb, si,
                 herm_hs=True, herm_c=False, m_less_n=True,
                 mtype=float, ham_ssq_q=False):
        """Initialization of the QuantumDot class."""
        self.si = si
        self.mtype = mtype
        self.ham_ssq_q = ham_ssq_q
        self.herm_hs = herm_hs
        self.herm_c = herm_c
        self.m_less_n = m_less_n
        self.hsingle = make_hsingle_dict(self, hsingle)
        self.coulomb = make_coulomb_dict(self, coulomb)
        self._init_hamiltonian()

    def _init_hamiltonian(self):
        si = self.si
        if si.indexing == 'sz':
            self.valslst = empty_szlst(si.nsingle, True)
            self.vecslst = empty_szlst(si.nsingle, True)
            self.hamlst = empty_szlst(si.nsingle, True)
        elif si.indexing == 'ssq':
            self.valslst = empty_ssqlst(si.nsingle, True)
            self.vecslst = empty_ssqlst(si.nsingle, True)
            self.hamlst = empty_szlst(si.nsingle, True)
        else:
            self.valslst = [None]*si.ncharge
            self.vecslst = [None]*si.ncharge
            self.hamlst = [None]*si.ncharge
        self.add(self.hsingle, self.coulomb, False)
        self.Ea = np.zeros(si.nmany, dtype=float)
        self.Ea_ext = None

    def add(self, hsingle=None, coulomb=None, updateq=True):
        """
        Adds a value to single particle Hamiltonian and Coulomb matrix elements
        and correspondingly redefines many-body Hamiltonian.

        Parameters
        ----------
        hsingle : dict
            Dictionary describing what single-particle Hamiltonian values to add.
            For example, hsingle[(state1, state2)] = value to add.
        coulomb : dict
            Dictionary describing what Coulomb matrix element values to add.
            For example, coulomb[(state1, state2, state3, state4)] = value to add.
        updateq : bool
            Specifies if the values of the single particle Hamiltonian
            and the Coulomb matrix elements will be updated.
            The many-body Hamiltonian will be updates in either case.
        """
        if updateq:
            hsingle = {} if hsingle is None else make_hsingle_dict(self, hsingle)
            coulomb = {} if coulomb is None else make_coulomb_dict(self, coulomb)
            for j0 in hsingle:
                if j0 in self.hsingle:
                    self.hsingle[j0] += hsingle[j0]
                else:
                    self.hsingle.update({j0: hsingle[j0]})
            for j0 in coulomb:
                if j0 in self.coulomb:
                    self.coulomb[j0] += coulomb[j0]
                else:
                    self.coulomb.update({j0: coulomb[j0]})
        #
        if self.si.indexing == 'sz':
            for charge in range(self.si.ncharge):
                for sz in szrange(charge, self.si.nsingle):
                    szlst = self.si.get_lst(charge, sz)
                    szind = sz_to_ind(sz, charge, self.si.nsingle)
                    self.hamlst[charge][szind] = construct_ham_coulomb(self, coulomb, szlst,
                                                                       self.hamlst[charge][szind])
                    self.hamlst[charge][szind] = construct_ham_hopping(self, hsingle, szlst,
                                                                       self.hamlst[charge][szind])
        elif self.si.indexing == 'ssq':
            for charge in range(self.si.ncharge):
                szlow = -(charge % 2)
                szrng = szrange(charge, self.si.nsingle) if self.ham_ssq_q else [szlow]
                for sz in szrng:
                    szlst = self.si.get_lst(charge, sz)
                    szind = sz_to_ind(sz, charge, self.si.nsingle)
                    self.hamlst[charge][szind] = construct_ham_coulomb(self, coulomb, szlst,
                                                                       self.hamlst[charge][szind])
                    self.hamlst[charge][szind] = construct_ham_hopping(self, hsingle, szlst,
                                                                       self.hamlst[charge][szind])
        else:
            for charge in range(self.si.ncharge):
                chargelst = self.si.get_lst(charge)
                self.hamlst[charge] = construct_ham_coulomb(self, coulomb, chargelst,
                                                            self.hamlst[charge])
                self.hamlst[charge] = construct_ham_hopping(self, hsingle, chargelst,
                                                            self.hamlst[charge])

    def change(self, hsingle=None, coulomb=None, updateq=True):
        """
        Changes the values of the single particle Hamiltonian and Coulomb matrix elements
        and correspondingly redefines many-body Hamiltonian.

        Parameters
        ----------
        hsingle : dict
            Dictionary describing what single-particle Hamiltonian values to add.
            For example, hsingle[(state1, state2)] = value to add.
        coulomb : dict
            Dictionary describing what Coulomb matrix element values to add.
            For example, coulomb[(state1, state2, state3, state4)] = value to add.
        updateq : bool
            Specifies if the values of the single particle Hamiltonian
            and the Coulomb matrix elements will be updated.
            The many-body Hamiltonian will be updates in either case.
        """
        hsingle = {} if hsingle is None else make_hsingle_dict(self, hsingle, True)
        coulomb = {} if coulomb is None else make_coulomb_dict(self, coulomb)
        #
        hsingle_add = {}
        coulomb_add = {}
        # Find the differences from the previous single-particle Hamiltonian
        # and from the previous Coulomb matrix elements
        for j0 in hsingle:
            if j0 in self.hsingle:
                hsingle_diff = hsingle[j0]-self.hsingle[j0]
                if hsingle_diff != 0:
                    hsingle_add.update({j0: hsingle_diff})
                    if updateq:
                        self.hsingle[j0] += hsingle_diff
            else:
                hsingle_diff = hsingle[j0]
                if hsingle_diff != 0:
                    hsingle_add.update({j0: hsingle_diff})
                    if updateq:
                        self.hsingle.update({j0: hsingle_diff})
        for j0 in coulomb:
            if j0 in self.coulomb:
                coulomb_diff = coulomb[j0]-self.coulomb[j0]
                if coulomb_diff != 0:
                    coulomb_add.update({j0: coulomb_diff})
                    if updateq:
                        self.coulomb[j0] += coulomb_diff
            else:
                coulomb_diff = coulomb[j0]
                if coulomb_diff != 0:
                    coulomb_add.update({j0: coulomb_diff})
                    if updateq:
                        self.coulomb.update({j0: coulomb_diff})
        # Then add the differences
        self.add(hsingle_add, coulomb_add, False)

    def diagonalise_charge(self, charge=None, sz=None, ssq=None):
        """Diagonalises Hamiltonian for given charge.
           Additionally, spin projection 'sz' and total spin 'ssq' can be specified
           if those symmetries exist for considered model.

        Parameters
        ----------
        charge : int
            Charge under consideration.
        sz : int
            Spin projection.
        ssq : int
            Total spin.
        """
        if charge is None:
            return None
        elif sz is None:
            chargelst = self.si.get_lst(charge)
            (self.valslst[charge],
             self.vecslst[charge]) = (
                construct_manybody_eigenstates(self, self.hsingle, self.coulomb,
                                               chargelst, self.hamlst[charge]))
            return self.valslst[charge], self.vecslst[charge]
        elif ssq is None and (self.si.indexing == 'sz' or self.si.indexing == 'ssq'):
            szlst = self.si.get_lst(charge, sz)
            szind = sz_to_ind(sz, charge, self.si.nsingle)
            (self.valslst[charge][szind],
             self.vecslst[charge][szind]) = (
                construct_manybody_eigenstates(self, self.hsingle, self.coulomb,
                                               szlst, self.hamlst[charge][szind]))
            return self.valslst[charge][szind], self.vecslst[charge][szind]
        elif self.si.indexing == 'ssq':
            szind = sz_to_ind(sz, charge, self.si.nsingle)
            ssqind = ssq_to_ind(ssq, sz)
            (self.valslst[charge][szind][ssqind],
             self.vecslst[charge][szind][ssqind]) = (
                construct_manybody_eigenstates_ssq(self, charge, sz, ssq,
                                                   self.hsingle, self.coulomb,
                                                   self.hamlst[charge][szind]))
            return self.valslst[charge][szind][ssqind], self.vecslst[charge][szind][ssqind]
        else:
            print("WARNING: No indexing by 'sz' or 'ssq'. Return None.")
            return None

    def diagonalise(self):
        """Diagonalises Hamiltonians for all charge states."""
        if self.si.indexing == 'sz':
            for charge in range(self.si.ncharge):
                for sz in szrange(charge, self.si.nsingle):
                    self.diagonalise_charge(charge, sz)
        elif self.si.indexing == 'ssq':
            for charge in range(self.si.ncharge):
                szlow = -(charge % 2)
                szind = sz_to_ind(szlow, charge, self.si.nsingle)
                (self.valslst[charge], self.vecslst[charge]) = (
                    construct_manybody_eigenstates_ssq_all(self, charge, self.hsingle, self.coulomb,
                                                           self.hamlst[charge][szind]))
        else:
            for charge in range(self.si.ncharge):
                self.diagonalise_charge(charge)
        self.set_Ea()

    def set_Ea(self):
        """Sets the many-body eigenstates using construct_Ea_manybody()."""
        self.Ea[:] = construct_Ea_manybody(self.valslst, self.si) class LeadsTunneling(object):
    """
    Class for defining tunneling amplitude from leads to the quantum dot.

    Attributes
    ----------
    nleads : int
        Number of the leads.
    tleads : dict, list or ndarray
        Dictionary, list or numpy array defining single particle tunneling amplitudes.
        numpy array has to be nleads by nsingle.
    si : StateIndexing
        StateIndexing or StateIndexingDM object.
    mulst : dict, list or ndarray
        Dictionary, list or numpy array containing chemical potentials of the leads.
    tlst : dict, list or ndarray
        Dictionary, list or numpy array containing temperatures of the leads.
    dlst : dict, list or ndarray
        Dictionary, list or numpy array containing bandwidths of the leads.
    mtype : type
        Defines type of Tba0 and Tba matrices. For example, float, complex, etc.
    Tba0 : ndarray
        nleads by nmany by nmany array, which contains many-body tunneling amplitude matrix
        in Fock basis.
    Tba : ndarray
        nleads by nmany by nmany array, which contains many-body tunneling amplitude matrix,
        which is used in calculations.
    """

    def __init__(self, nleads, tleads, si, mulst, tlst, dlst, mtype=complex):
        """Initialization of the LeadsTunneling class."""
        si.nleads = nleads
        si.nleads_sym = nleads//2 if si.symmetry == 'spin' else nleads
        #
        self.si = si
        self.tleads = make_tleads_dict(tleads, si)
        self.tleads_array = make_tleads_array(self.tleads, si)
        self.mulst = make_array(None, mulst, si)
        self.tlst = make_array(None, tlst, si)
        self.dlst = make_array_dlst(None, dlst, si)
        self.mtype = mtype
        self._init_coupling()

    def _init_coupling(self):
        self.Tba0 = construct_Tba(self, self.tleads)
        self.Tba = np.array(self.Tba0)

    def add(self, tleads=None, mulst=None, tlst=None, dlst=None, updateq=True, lstq=True):
        """
        Adds a value to single particle tunneling amplitudes and correspondingly redefines
        many-body tunneling matrix Tba.

        Parameters
        ----------
        tleads,mulst,tlst,dlst : dict, list, int
            Dictionaries describing what values to add.
            For example, tleads[(lead, state)] = value to add.
        updateq : bool
            Specifies if the values of the single particle amplitudes will be updated.
            The many-body tunneling amplitudes Tba will be updates in either case.
        lstq : bool
            Determines if the values will be added to mulst, tlst, dlst.
        """
        if lstq:
            if mulst is not None:
                self.mulst += make_array(None, mulst, self.si)
            if tlst is not None:
                self.tlst += make_array(None, tlst, self.si)
            if dlst is not None:
                self.dlst += make_array_dlst(None, dlst, self.si)
        if tleads is not None:
            if updateq:
                tleads = make_tleads_dict(tleads, self.si)
                for j0 in tleads:
                    if j0 in self.tleads:
                        self.tleads[j0] += tleads[j0]
                    else:
                        self.tleads.update({j0: tleads[j0]})
                self.tleads_array = make_tleads_array(self.tleads, self.si)
            self.Tba0 = construct_Tba(self, tleads, self.Tba0)

    def change(self, tleads=None, mulst=None, tlst=None, dlst=None):
        """
        Changes the values of the single particle tunneling amplitudes and correspondingly redefines
        many-body tunneling matrix Tba.

        Parameters
        ----------
        tleads,mulst,tlst,dlst : dict, list, int
            Dictionaries describing what values to change.
            For example, tleads[(lead, state)] = the new value.
        """
        if mulst is not None:
            self.mulst[:] = make_array(self.mulst, mulst, self.si)
        if tlst is not None:
            self.tlst[:] = make_array(self.tlst, tlst, self.si)
        if dlst is not None:
            self.dlst[:] = make_array_dlst(self.dlst, dlst, self.si)
        #
        if tleads is not None:
            tleads = make_tleads_dict(tleads, self.si, True)
            # Find the differences from the previous tunneling amplitudes
            tleads_add = {}
            for j0 in tleads:
                if j0 in self.tleads:
                    tleads_diff = tleads[j0]-self.tleads[j0]
                    if tleads_diff != 0:
                        tleads_add.update({j0: tleads_diff})
                        self.tleads[j0] += tleads_diff
                else:
                    tleads_diff = tleads[j0]
                    if tleads_diff != 0:
                        tleads_add.update({j0: tleads_diff})
                        self.tleads.update({j0: tleads_diff})
            # Add the differences
            self.add(tleads_add, updateq=False, lstq=False)

    def rotate(self, vecslst, indexing=None):
        """
        Rotates tunneling amplitude matrix Tba0 in Fock basis to Tba,
        which is in eigenstate basis of the quantum dot.

        Parameters
        ----------
        vecslst : list of ndarrays
            List of size ncharge containing arrays defining eigenvector matrices for given charge.
        indexing : string
            Specifies what kind of rotation procedure to use. Default is si.indexing.
        """
        self.Tba[:] = rotate_Tba(self.Tba0, vecslst, self.si, indexing, self.mtype)

    def use_Tba0(self):
        """
        Sets the Tba matrix for calculation to Tba0 in the Fock basis.
        """
        self.Tba[:] = self.Tba0

    def update_Tba0(self, nleads, tleads, mtype=complex):
        """
        Updates the Tba0 in the Fock basis using new single-particle tunneling amplitudes.

        Parameters
        ----------
        nleads : int
            Number of leads.
        tleads : dict
            The new single-particle tunneling amplitudes. See attribute tleads.
        mtype : type
            Defines type of Tba0 and Tba matrices. For example, float, complex, etc.
        """
        si = self.si
        si.nleads = nleads
        si.nleads_sym = nleads//2 if si.symmetry == 'spin' else nleads
        self.tleads = make_tleads_dict(tleads, si)
        self.tleads_array = make_tleads_array(tleads, self.si)
        self.mtype = mtype
        self.Tba0 = construct_Tba(self, tleads) class FunctionProperties(object):
    """
    Class containing miscellaneous variables for Approach and Approach2vN classes.

    Attributes
    ----------
    symq : bool
        For symq=False keep all equations in the kernel, and the matrix is of size N by N+1.
        For symq=True replace one equation by the normalisation condition,
        and the matrix is square N by N.
    norm_row : int
        If symq=True this row will be replaced by normalisation condition in the kernel matrix.
    solmethod : string
        String specifying the solution method of the equation L(Phi0)=0.
        The possible values are matrix inversion 'solve' and least squares 'lsqr'.
        Method 'solve' works only when symq=True.
        For matrix free methods (used when mfreeq=True) the possible values are
        'krylov', 'broyden', etc.
    itype : int
        Type of integral for first order approach calculations.
        itype=0: the principal parts are evaluated using Fortran integration package QUADPACK \
                 routine dqawc through SciPy.
        itype=1: the principal parts are kept, but approximated by digamma function valid for \
                 large bandwidth D.
        itype=2: the principal parts are neglected.
        itype=3: the principal parts are neglected and infinite bandwidth D is assumed.
    dqawc_limit : int
        For itype=0 dqawc_limit determines the maximum number of sub-intervals
        in the partition of the given integration interval.
    mfreeq : bool
        If mfreeq=True the matrix free solution method is used for first order methods.
    phi0_init : array
        For mfreeq=True the initial value of zeroth order density matrix elements.
    mtype_qd : float or complex
        Type for the many-body quantum dot Hamiltonian matrix.
    mtype_leads : float or complex
        Type for the many-body tunneling matrix Tba.
    kpnt_left, kpnt_right : int
        Number of points Ek_grid is extended to the left and the right for '2vN' approach.
    ht_ker : array
        Kernel used when performing Hilbert transform using FFT.
        It is generated using specfunc.kernel_fredriksen(n).
    emin, emax : float
        Minimal and maximal energy in the updated Ek_grid generated by neumann2py.get_grid_ext(sys).
        Note that emin<=Dmin and emax>=Dmax.
    dmin, dmax : float
        Bandedge Dmin and Dmax values of the lead electrons.
    ext_fct : float
        Multiplication factor used in neumann2py.get_grid_ext(sys), when determining emin and emax.
    suppress_err : bool
        Determines whether to print the warning when the inversion of the kernel failed.
    off_diag_corrections: bool
        Determines wether to include first oder off-diagonal corrections to the kernel in the
        RTD approach.
    """

    def __init__(self,
                 kerntype='2vN', symq=True, norm_row=0, solmethod=None,
                 itype=0, dqawc_limit=10000, mfreeq=False, phi0_init=None,
                 mtype_qd=float, mtype_leads=complex, kpnt=None, dband=None,
                 off_diag_corrections=True):
        self.kerntype = kerntype
        self.symq = symq
        self.norm_row = norm_row
        self.solmethod = solmethod
        #
        self.itype = itype
        self.dqawc_limit = dqawc_limit
        #
        self.mfreeq = mfreeq
        self.phi0_init = phi0_init
        #
        self.mtype_qd = mtype_qd
        self.mtype_leads = mtype_leads
        #
        self.kpnt = kpnt
        self.dband = dband
        #
        self.kpnt_left = 0
        self.kpnt_right = 0
        self.ht_ker = None
        #
        self.dmin, self.dmax = 0, 0
        self.emin, self.emax = 0, 0
        self.ext_fct = 1.1
        #
        self.suppress_err = False
        self.suppress_wrn = [False]
        #
        self.off_diag_corrections = off_diag_corrections

    def print_error(self, exept):
        if not self.suppress_err:
            print("WARNING: Could not solve the linear set of equations.\n" +
                  "  Error from the solver: " + str(exept) + "\n"
                  "  The reasons for such a failure can be various:\n" +
                  "  1. Some of the transport channels may be outside the bandwidth D of the leads.\n" +
                  "     In this case removing some of the states with the method [remove_states()] will help.\n" +
                  "  2. Replacement of one of the equations with the normalisation condition.\n" +
                  "     In this case try to use different [norm_row]\n"+
                  "     or solve the linear system using [symq=False] and the solution method [solmethod='lsqr'].\n"
                  "  This warning will not be shown again.\n"
                  "  To check if the solution succeeded check the property [success].")
            self.suppress_err = True

    def print_warning(self, i, message):
        if not self.suppress_wrn[i]:
            print(message)
            self.suppress_wrn[i] = True def get_phi0(self, b_, bp_):
    """
    Get the reduced density matrix element corresponding to
    many-body states b and bp.

    Parameters
    ----------
    self : Builder or Approach
        The system given as Builder or Approach object.
    b_,bp_ : int
        Labels of the many-body states.

    Returns
    --------
    phi0bbp : complex
        A matrix element of the reduced density matrix (complex number).
    """
    b = self.si.states_order[b_]
    bp = self.si.states_order[bp_]
    bcharge = sum(self.si.get_state(b))
    bpcharge = sum(self.si.get_state(bp))
    phi0bbp = 0.0
    if self.funcp.kerntype == 'Pauli':
        if b == bp:
            ind = self.si.get_ind_dm0(b, b, bcharge, maptype=1)
            phi0bbp = self.phi0[ind]
    elif bcharge == bpcharge:
        ind = self.si.get_ind_dm0(b, bp, bcharge, maptype=1)
        conj = self.si.get_ind_dm0(b, bp, bcharge, maptype=3)
        if ind != -1:
            if type(self.si).__name__ == 'StateIndexingDMc':
                phi0bbp = self.phi0[ind]
            else:
                ndm0, npauli = self.si.ndm0, self.si.npauli
                phi0bbp = (self.phi0[ind] + 1j*self.phi0[ndm0-npauli+ind]
                           * (+1 if conj else -1)
                           * (0 if ind < npauli else 1))
    return phi0bbp def get_phi1(self, l, c_, b_):
    """
    Get the energy integrated current amplitudes corresponding to
    lead l and many-body states c and b.

    Parameters
    ----------
    self : Builder or Approach
        The system given as Builder or Approach object.
    l : int
        Label of the lead channel.
    c_,b_ : int
        Labels of the many-body states.

    Returns
    --------
    phi0bbp : complex
        A matrix element of the reduced density matrix (complex number).
    """
    if self.funcp.kerntype == 'Pauli':
        return None
    else:
        c = self.si.states_order[c_]
        b = self.si.states_order[b_]
        ccharge = sum(self.si.get_state(c))
        bcharge = sum(self.si.get_state(b))
        phi1cb = 0.0
        if ccharge == bcharge+1:
            ind = self.si.get_ind_dm1(c, b, bcharge)
            phi1cb = self.phi1[l, ind]
        elif ccharge+1 == bcharge:
            ind = self.si.get_ind_dm1(b, c, ccharge)
            phi1cb = self.phi1[l, ind].conjugate()
        return phi1cb def print_state(self, b_, eps=0.0, prntq=True, filename=None, separator=''):
    """
    Prints properties of given many-body eigenstate of the quantum dot Hamiltonian

    Parameters
    ----------
    self : Builder or Approach
        The system given as Builder or Approach object.
    b_ : int
        Index of the many-body eigenstate.
    eps : float
        Value at which expansion coefficients of the many-body eigenstate in Fock basis
        are neglected.
    prntq : bool
        If true then eigenstate properties are printed into command line.
    filename : string
        File to which to print eigenstate properties. The output is appended to this file.
    separator : string
        String added at the beginning of the output.
    """
    si = self.si
    b = si.states_order[b_]
    Eb = self.qd.Ea[b]
    sout = ('Energy: ' + str(Eb)+'\n' +
            'Eigenstate coefficients:')
    #
    if si.indexing == 'sz':
        charge, sz, alpha = si.ind_qn[b]
        sout = 'S^z [in units of 1/2]: ' + str(sz) + '\n' + sout
        szind = sz_to_ind(sz, charge, si.nsingle)
        ind = si.szlst[charge][szind].index(b)
        coeffs = np.array([abs(self.qd.vecslst[charge][szind][:, ind]),
                           range(len(si.szlst[charge][szind]))])
        coeffs = np.array(multiarray_sort(coeffs)[1], dtype=int)
        # for j1 in range(len(si.szlst[charge][szind])):
        for j1 in reversed(coeffs):
            sn = si.szlst[charge][szind][j1]
            val = self.qd.vecslst[charge][szind][j1, ind]
            if abs(val) >= eps:
                sout = sout + '\n'+'    |'+si.get_state(sn, strq=True)+'>: ' + str(val)
    elif si.indexing == 'ssq':
        charge, sz, ssq, alpha = si.ind_qn[b]
        sout = ('S^z [in units of 1/2]: ' + str(sz) + '\n' +
                'S^2 [in units of 1/2]: ' + str(ssq) + '\n' + sout)
        szind = sz_to_ind(sz, charge, si.nsingle)
        ssqind = ssq_to_ind(ssq, sz)
        ind = si.ssqlst[charge][szind][ssqind].index(b)
        coeffs = np.array([abs(self.qd.vecslst[charge][szind][ssqind][:, ind]),
                           range(len(si.szlst[charge][szind]))])
        coeffs = np.array(multiarray_sort(coeffs)[1], dtype=int)
        # for j1 in range(len(si.szlst[charge][szind])):
        for j1 in reversed(coeffs):
            sn = si.szlst[charge][szind][j1]
            val = self.qd.vecslst[charge][szind][ssqind][j1, ind]
            if abs(val) >= eps:
                sout = sout + '\n'+'    |'+si.get_state(sn, strq=True)+'>: ' + str(val)
    else:
        charge, alpha = si.ind_qn[b]
        # ind = si.dictdm[b]
        ind = si.chargelst[charge].index(b)
        coeffs = np.array([abs(self.qd.vecslst[charge][:, ind]),
                           range(len(si.chargelst[charge]))])
        coeffs = np.array(multiarray_sort(coeffs)[1], dtype=int)
        # for j1 in range(len(si.chargelst[charge])):
        for j1 in reversed(coeffs):
            sn = si.chargelst[charge][j1]
            val = self.qd.vecslst[charge][j1, ind]
            if abs(val) >= eps:
                sout = sout + '\n'+'    |'+si.get_state(sn, strq=True)+'>: ' + str(val)
    #
    sout = (separator +
            'Sorted label: ' + str(b_) + '\n' +
            'Original label: ' + str(b) + '\n' +
            'Charge: ' + str(charge) + '\n' + sout)
    #
    if filename is not None:
        with open(filename, 'a') as f:
            f.writelines(sout+'\n')
    if prntq:
        print(sout)
    pass def print_all_states(self, filename, eps=0.0, separator='', mode='w'):
    """
    Prints properties of all many-body eigenstates to a file.

    Parameters
    ----------
    self : Builder or Approach
        The system given as Builder or Approach object.
    filename : string
        File to which to print eigenstate properties. The output is appended to this file.
    eps : float
        Value at which expansion coefficients of the many-body eigenstate in Fock basis
        are neglected.
    separator : string
        String added at the beginning of the output.
    mode : string
        Mode in which the file will be treated. For example, to append use 'a' or
        to write into clea file use 'w'.
    """
    with open(filename, mode) as f:
        f.close()
    for ind in range(self.si.nmany):
        print_state(self, ind, eps, False, filename, separator)
    pass def sort_eigenstates(self, srt=None):
    """
    Sort many-body states of the system by given order of properties.

    Parameters
    ----------
    self : Builder or Approach.
        The system given as Approach, Approach2vN, or Builder object.
    srt : list
        List specifying in which order of properties to sort.
        For example, in the case of 'ssq' indexing we have  such convention:
        0 - energy
        1 - charge
        2 - spin projection :math:`S_{z}`
        3 - total spin :math:`S^{2}`
        The default sorting order for 'ssq' indexing is srt=[1, 2, 3, 0]

    self.qd.Ea_ext : ndarray
        (Modifies) A two-dimensional numpy array containing in the zeroth row energies,
        first row charge, and etc.
    self.si.states_order : ndarray
        (Modifies) A numpy row containing a new order of many-body states.
    """
    if srt is not None:
        if self.qd.Ea_ext is None:
            self.qd.Ea_ext = construct_Ea_extended(self)
        else:
            self.qd.Ea_ext[0] = self.qd.Ea
        srt.append(-1)
        self.si.states_order = np.array(multiarray_sort(self.qd.Ea_ext, srt)[-1], dtype=int)
    else:
        self.si.states_order = range(self.si.nmany) def remove_coherences(self, dE):
    """
    Remove the coherences with energy difference larger than dE.

    Parameters
    ----------
    self : Builder or Approach
        The system given as Builder or Approach.
    dE : float
        Energy difference.

    self.si.mapdm0 : list
        (Modifies) List showing which density matrix elements are mapped to each other due to symmetries
        and which density matrix elements are neglected (entries with values -1).
    """
    # Find which coherences to remove
    si, E = self.si, self.qd.Ea
    ilst, indlst = [], []
    for i in range(si.ndm0_):
        if si.mapdm0[i] != -1:
            ind = si.mapdm0[i]
            b, bp = si.inddm0[ind]
            if abs(E[b]-E[bp]) > dE:
                si.mapdm0[i] = -1
            else:
                ilst.append(i)
                indlst.append(ind)
    # Count the new number of density matrix elements
    srt = multiarray_sort(np.array([ilst, indlst]), [1, 0])
    ilst, indlst = srt[0], srt[1]
    ind, count = 0, 0
    for i in range(len(indlst)):
        if indlst[i] == ind:
            indlst[i] = count
        else:
            count += 1
            ind = indlst[i]
            indlst[i] = count
    # Relabel the density matrix elements
    si.ndm0 = count+1
    si.ndm0r = 2*si.ndm0-si.npauli
    for i in range(len(indlst)):
        si.mapdm0[ilst[i]] = indlst[i] def remove_states(self, dE):
    """
    Remove the states with energy dE larger than the ground state
    for the transport calculations.

    Parameters
    ----------
    self : Builder or Approach
        The system given as Builder or Approach object.
    dE : float
        Energy above the ground state.

    self.si.statesdm : list
        (Modifies) List containing indices of many-body state under consideration.
    """
    Emax = min(self.qd.Ea)+dE
    statesdm = [[] for _ in range(self.si.nsingle+1)]
    for charge in range(self.si.ncharge):
        for b in self.si.chargelst[charge]:
            if self.qd.Ea[b] < Emax:
                statesdm[charge].append(b)
    self.si.set_statesdm(statesdm) def use_all_states(self):
    """
    Use all states for the transport calculations.

    Parameters
    ----------
    self : Builder or Approach
        The system given as Builder or Approach.

    self.si.statesdm : list
        (Modifies) List containing indices of many-body state under consideration.
    """
    self.si.set_statesdm(self.si.chargelst) def validate_kerntype(kerntype):
    if isinstance(kerntype, str):
        if kerntype not in {'Pauli', 'Lindblad', 'Redfield', '1vN', '2vN', 'pyPauli',
                    'pyLindblad', 'pyRedfield', 'py1vN', 'py2vN', 'pyRTD', 'RTD'}:
            print("WARNING: Allowed kerntype values are: " +
                  "\'Pauli\', \'Lindblad\', \'Redfield\', \'1vN\', \'2vN\', " +
                  "\'pyPauli\', \'pyLindblad\', \'pyRedfield\', \'py1vN\', \'py2vN\', \'RTD\'. " +
                  "Using default kerntype=\'Pauli\'.")
            kerntype = 'Pauli'
    return kerntype def validate_itype(itype, kerntype):
    if kerntype not in ('RTD', 'pyRTD'):
        if itype not in {0, 1, 2, 3}:
            print("WARNING: itype needs to be 0, 1, 2, or 3. Using default itype=0.")
            itype = 0
    else:
        if itype != 1:
            print("WARNING: only itype=1 is supported by the RTD approach. Using itype=1.")
            itype = 1
    return itype def validate_indexing(indexing, symmetry, kerntype):
    if indexing is None:
        if symmetry == 'spin' and kerntype in {'pyRTD', 'RTD'}:
            print("WARNING: symmetry=\'spin\' is not supported by the RTD approach. " +
                  "Using default indexing=\'charge\'.")
            indexing = 'charge'
            symmetry = None
        elif symmetry == 'spin' and kerntype not in {'py2vN', '2vN'}:
            indexing = 'ssq'
        else:
            indexing = 'charge'

    if indexing not in {'Lin', 'charge', 'sz', 'ssq'}:
        print("WARNING: Allowed indexing values are: \'Lin\', \'charge\', \'sz\', \'ssq\'. " +
              "Using default indexing=\'charge\'.")
        indexing = 'charge'

    if indexing not in {'Lin', 'charge'} and kerntype in {'py2vN', '2vN'}:
        print("WARNING: For the 2vN approach indexing needs to be \'Lin\' or \'charge\'. " +
              "Using indexing=\'charge\' as a default.")
        indexing = 'charge'

    if indexing != 'charge' and kerntype in {'pyRTD', 'RTD'}:
        print("WARNING: For the RTD approach indexing needs to be \'charge\'. " +
              "Using indexing=\'charge\' as a default.")
        indexing = 'charge'

    return indexing, symmetry class ApproachPauli(Approach):

    kerntype = 'pyPauli'

    def get_kern_size(self):
        return self.si.npauli

    def prepare_arrays(self):
        Approach.prepare_arrays(self)
        nleads, ndm1 = self.si.nleads, self.si.ndm1
        self.paulifct = np.zeros((nleads, ndm1, 2), dtype=doublenp)

    def clean_arrays(self):
        Approach.clean_arrays(self)
        self.paulifct.fill(0.0)

    def generate_fct(self):
        """
        Make factors used for generating Pauli master equation kernel.

        Parameters
        ----------
        paulifct : array
            (Modifies) Factors used for generating Pauli master equation kernel.
        """
        E, Tba, si = self.qd.Ea, self.leads.Tba, self.si
        mulst, tlst, dlst = self.leads.mulst, self.leads.tlst, self.leads.dlst
        ncharge, nleads, statesdm = si.ncharge, si.nleads, si.statesdm

        itype = self.funcp.itype
        paulifct = self.paulifct
        for charge in range(ncharge-1):
            ccharge = charge+1
            bcharge = charge
            for c, b in itertools.product(statesdm[ccharge], statesdm[bcharge]):
                cb = si.get_ind_dm1(c, b, bcharge)
                Ecb = E[c]-E[b]
                for l in range(nleads):
                    xcb = (Tba[l, b, c]*Tba[l, c, b]).real
                    rez = func_pauli(Ecb, mulst[l], tlst[l], dlst[l, 0], dlst[l, 1], itype)
                    paulifct[l, cb, 0] = xcb*rez[0]
                    paulifct[l, cb, 1] = xcb*rez[1]

    def generate_kern(self):
        """
        Generate Pauli master equation kernel.

        Parameters
        ----------
        kern : array
            (Modifies) Kernel matrix for Pauli master equation.
        """
        si, kh = self.si, self.kernel_handler
        ncharge, statesdm = si.ncharge, si.statesdm

        for bcharge in range(ncharge):
            for b in statesdm[bcharge]:
                if not kh.is_unique(b, b, bcharge):
                    continue
                self.generate_coupling_terms(b, b, bcharge)

    def generate_coupling_terms(self, b, bp, bcharge):
        paulifct = self.paulifct
        si, kh = self.si, self.kernel_handler
        nleads, statesdm = si.nleads, si.statesdm

        acharge = bcharge-1
        ccharge = bcharge+1

        bb = si.get_ind_dm0(b, b, bcharge)
        for a in statesdm[acharge]:
            aa = si.get_ind_dm0(a, a, acharge)
            ba = si.get_ind_dm1(b, a, acharge)
            fctm, fctp = 0, 0
            for l in range(nleads):
                fctm -= paulifct[l, ba, 1]
                fctp += paulifct[l, ba, 0]
            kh.set_matrix_element_pauli(fctm, fctp, bb, aa)
        for c in statesdm[ccharge]:
            cc = si.get_ind_dm0(c, c, ccharge)
            cb = si.get_ind_dm1(c, b, bcharge)
            fctm, fctp = 0, 0
            for l in range(nleads):
                fctm -= paulifct[l, cb, 0]
                fctp += paulifct[l, cb, 1]
            kh.set_matrix_element_pauli(fctm, fctp, bb, cc)

    def generate_current(self):
        """
        Calculates currents using Pauli master equation approach.

        Parameters
        ----------
        current : array
            (Modifies) Values of the current having nleads entries.
        energy_current : array
            (Modifies) Values of the energy current having nleads entries.
        heat_current : array
            (Modifies) Values of the heat current having nleads entries.
        """
        phi0, E, paulifct, si = self.phi0, self.qd.Ea, self.paulifct, self.si
        ncharge, nleads, statesdm = si.ncharge, si.nleads, si.statesdm

        current = self.current
        energy_current = self.energy_current

        for charge in range(ncharge-1):
            ccharge = charge+1
            bcharge = charge
            for c in statesdm[ccharge]:
                cc = si.get_ind_dm0(c, c, ccharge)
                for b in statesdm[bcharge]:
                    bb = si.get_ind_dm0(b, b, bcharge)
                    cb = si.get_ind_dm1(c, b, bcharge)
                    for l in range(nleads):
                        fct1 = +phi0[bb]*paulifct[l, cb, 0]
                        fct2 = -phi0[cc]*paulifct[l, cb, 1]
                        current[l] += fct1 + fct2
                        energy_current[l] += -(E[b]-E[c])*(fct1 + fct2)

        self.heat_current[:] = energy_current - current*self.leads.mulst class ApproachLindblad(Approach):

    kerntype = 'pyLindblad'

    def prepare_arrays(self):
        Approach.prepare_arrays(self)
        Tba, mtype = self.leads.Tba, self.leads.mtype
        self.tLba = np.zeros(Tba.shape, dtype=mtype)

    def clean_arrays(self):
        Approach.clean_arrays(self)
        self.tLba.fill(0.0)

    def generate_fct(self):
        """
        Make factors used for generating Lindblad master equation kernel.

        Parameters
        ----------
        tLba : array
            (Modifies) Jump operator matrix in many-body basis.
        """
        Tba, E, si = self.leads.Tba, self.qd.Ea, self.si
        mulst, tlst, dlst = self.leads.mulst, self.leads.tlst, self.leads.dlst
        itype = self.funcp.itype
        ncharge, nleads, statesdm = si.ncharge, si.nleads, si.statesdm

        tLba = self.tLba
        for charge in range(ncharge-1):
            bcharge = charge+1
            acharge = charge
            for b, a in itertools.product(statesdm[bcharge], statesdm[acharge]):
                Eba = E[b]-E[a]
                for l in range(nleads):
                    fct1, fct2 = func_pauli(Eba, mulst[l], tlst[l], dlst[l, 0], dlst[l, 1], itype)
                    tLba[l, b, a] = np.sqrt(fct1)*Tba[l, b, a]
                    tLba[l, a, b] = np.sqrt(fct2)*Tba[l, a, b]

    def generate_coupling_terms(self, b, bp, bcharge):
        tLba = self.tLba
        si, kh = self.si, self.kernel_handler
        nleads, statesdm = si.nleads, si.statesdm

        acharge = bcharge-1
        ccharge = bcharge+1

        for a, ap in itertools.product(statesdm[acharge], statesdm[acharge]):
            if kh.is_included(a, ap, acharge):
                fct_aap = 0
                for l in range(nleads):
                    fct_aap += tLba[l, b, a]*tLba[l, bp, ap].conjugate()
                kh.set_matrix_element(1j*fct_aap, b, bp, bcharge, a, ap, acharge)
        # --------------------------------------------------
        for bpp in statesdm[bcharge]:
            if kh.is_included(bpp, bp, bcharge):
                fct_bppbp = 0
                for a in statesdm[acharge]:
                    for l in range(nleads):
                        fct_bppbp += -0.5*tLba[l, a, b].conjugate()*tLba[l, a, bpp]
                for c in statesdm[ccharge]:
                    for l in range(nleads):
                        fct_bppbp += -0.5*tLba[l, c, b].conjugate()*tLba[l, c, bpp]
                kh.set_matrix_element(1j*fct_bppbp, b, bp, bcharge, bpp, bp, bcharge)
            # --------------------------------------------------
            if kh.is_included(b, bpp, bcharge):
                fct_bbpp = 0
                for a in statesdm[acharge]:
                    for l in range(nleads):
                        fct_bbpp += -0.5*tLba[l, a, bpp].conjugate()*tLba[l, a, bp]
                for c in statesdm[ccharge]:
                    for l in range(nleads):
                        fct_bbpp += -0.5*tLba[l, c, bpp].conjugate()*tLba[l, c, bp]
                kh.set_matrix_element(1j*fct_bbpp, b, bp, bcharge, b, bpp, bcharge)
        # --------------------------------------------------
        for c, cp in itertools.product(statesdm[ccharge], statesdm[ccharge]):
            if kh.is_included(c, cp, ccharge):
                fct_ccp = 0
                for l in range(nleads):
                    fct_ccp += tLba[l, b, c]*tLba[l, bp, cp].conjugate()
                kh.set_matrix_element(1j*fct_ccp, b, bp, bcharge, c, cp, ccharge)
        # --------------------------------------------------

    def generate_current(self):
        """
        Calculates currents using Lindblad approach.

        Parameters
        ----------
        current : array
            (Modifies) Values of the current having nleads entries.
        energy_current : array
            (Modifies) Values of the energy current having nleads entries.
        heat_current : array
            (Modifies) Values of the heat current having nleads entries.
        """
        phi0p, E, tLba, si = self.phi0, self.qd.Ea, self.tLba, self.si
        ndm0, npauli = si.ndm0, si.npauli
        ncharge, nleads, statesdm = si.ncharge, si.nleads, si.statesdm

        current = self.current
        energy_current = self.energy_current

        kh = self.kernel_handler
        for charge in range(ncharge):
            ccharge = charge+1
            bcharge = charge
            acharge = charge-1

            for b, bp in itertools.product(statesdm[bcharge], statesdm[bcharge]):
                if not kh.is_included(b, bp, bcharge):
                    continue
                phi0bbp = kh.get_phi0_element(b, bp, bcharge)

                for l in range(nleads):
                    current_l, energy_current_l = 0, 0

                    for a in statesdm[acharge]:
                        fcta = tLba[l, a, b]*phi0bbp*tLba[l, a, bp].conjugate()
                        current_l -= fcta
                        energy_current_l += (E[a]-0.5*(E[b]+E[bp]))*fcta
                    for c in statesdm[ccharge]:
                        fctc = tLba[l, c, b]*phi0bbp*tLba[l, c, bp].conjugate()
                        current_l += fctc
                        energy_current_l += (E[c]-0.5*(E[b]+E[bp]))*fctc

                    current[l] += current_l.real
                    energy_current[l] += energy_current_l.real

        self.heat_current[:] = energy_current - current*self.leads.mulst class ApproachRedfield(Approach):

    kerntype = 'pyRedfield'

    def prepare_arrays(self):
        Approach1vN.prepare_arrays(self)

    def clean_arrays(self):
        Approach1vN.clean_arrays(self)

    def generate_fct(self):
        Approach1vN.generate_fct(self)

    def generate_coupling_terms(self, b, bp, bcharge):
        Tba, phi1fct = self.leads.Tba, self.phi1fct
        si, kh = self.si, self.kernel_handler
        nleads, statesdm = si.nleads, si.statesdm

        acharge = bcharge-1
        ccharge = bcharge+1

        # --------------------------------------------------
        for a, ap in itertools.product(statesdm[acharge], statesdm[acharge]):
            if kh.is_included(a, ap, acharge):
                bpap = si.get_ind_dm1(bp, ap, acharge)
                ba = si.get_ind_dm1(b, a, acharge)
                fct_aap = 0
                for l in range(nleads):
                    fct_aap += (+ Tba[l, b, a]*Tba[l, ap, bp]*phi1fct[l, bpap, 0].conjugate()
                                - Tba[l, b, a]*Tba[l, ap, bp]*phi1fct[l, ba, 0])
                kh.set_matrix_element(fct_aap, b, bp, bcharge, a, ap, acharge)
        # --------------------------------------------------
        for bpp in statesdm[bcharge]:
            if kh.is_included(bpp, bp, bcharge):
                fct_bppbp = 0
                for a in statesdm[acharge]:
                    bppa = si.get_ind_dm1(bpp, a, acharge)
                    for l in range(nleads):
                        fct_bppbp += +Tba[l, b, a]*Tba[l, a, bpp]*phi1fct[l, bppa, 1].conjugate()
                for c in statesdm[ccharge]:
                    cbpp = si.get_ind_dm1(c, bpp, bcharge)
                    for l in range(nleads):
                        fct_bppbp += +Tba[l, b, c]*Tba[l, c, bpp]*phi1fct[l, cbpp, 0]
                kh.set_matrix_element(fct_bppbp, b, bp, bcharge, bpp, bp, bcharge)
            # --------------------------------------------------
            if kh.is_included(b, bpp, bcharge):
                fct_bbpp = 0
                for a in statesdm[acharge]:
                    bppa = si.get_ind_dm1(bpp, a, acharge)
                    for l in range(nleads):
                        fct_bbpp += -Tba[l, bpp, a]*Tba[l, a, bp]*phi1fct[l, bppa, 1]
                for c in statesdm[ccharge]:
                    cbpp = si.get_ind_dm1(c, bpp, bcharge)
                    for l in range(nleads):
                        fct_bbpp += -Tba[l, bpp, c]*Tba[l, c, bp]*phi1fct[l, cbpp, 0].conjugate()
                kh.set_matrix_element(fct_bbpp, b, bp, bcharge, b, bpp, bcharge)
        # --------------------------------------------------
        for c, cp in itertools.product(statesdm[ccharge], statesdm[ccharge]):
            if kh.is_included(c, cp, ccharge):
                cpbp = si.get_ind_dm1(cp, bp, bcharge)
                cb = si.get_ind_dm1(c, b, bcharge)
                fct_ccp = 0
                for l in range(nleads):
                    fct_ccp += (+ Tba[l, b, c]*Tba[l, cp, bp]*phi1fct[l, cpbp, 1]
                                - Tba[l, b, c]*Tba[l, cp, bp]*phi1fct[l, cb, 1].conjugate())
                kh.set_matrix_element(fct_ccp, b, bp, bcharge, c, cp, ccharge)
        # --------------------------------------------------

    def generate_current(self):
        E, Tba = self.qd.Ea, self.leads.Tba
        phi1fct, phi1fct_energy = self.phi1fct, self.phi1fct_energy

        si = self.si
        ncharge, nleads, statesdm = si.ncharge, si.nleads, si.statesdm

        phi1 = self.phi1
        current = self.current
        energy_current = self.energy_current

        kh = self.kernel_handler
        for charge in range(ncharge-1):
            ccharge = charge+1
            bcharge = charge
            for c, b in itertools.product(statesdm[ccharge], statesdm[bcharge]):
                cb = si.get_ind_dm1(c, b, bcharge)

                for l in range(nleads):
                    current_l, energy_current_l = 0, 0

                    for bp in statesdm[bcharge]:
                        if not kh.is_included(bp, b, bcharge):
                            continue
                        phi0bpb = kh.get_phi0_element(bp, b, bcharge)

                        cbp = si.get_ind_dm1(c, bp, bcharge)
                        fct1 = phi1fct[l, cbp, 0]
                        fct1h = phi1fct_energy[l, cbp, 0]

                        phi1[l, cb] += Tba[l, c, bp]*phi0bpb*fct1
                        current_l += Tba[l, b, c]*Tba[l, c, bp]*phi0bpb*fct1
                        energy_current_l += Tba[l, b, c]*Tba[l, c, bp]*phi0bpb*fct1h

                    for cp in statesdm[ccharge]:
                        if not kh.is_included(c, cp, ccharge):
                            continue
                        phi0ccp = kh.get_phi0_element(c, cp, ccharge)

                        cpb = si.get_ind_dm1(cp, b, bcharge)
                        fct2 = phi1fct[l, cpb, 1]
                        fct2h = phi1fct_energy[l, cpb, 1]

                        phi1[l, cb] += Tba[l, cp, b]*phi0ccp*fct2
                        current_l += Tba[l, b, c]*phi0ccp*Tba[l, cp, b]*fct2
                        energy_current_l += Tba[l, b, c]*phi0ccp*Tba[l, cp, b]*fct2h

                    current[l] += -2*current_l.imag
                    energy_current[l] += -2*energy_current_l.imag

        self.heat_current[:] = energy_current - current*self.leads.mulst class Approach1vN(Approach):

    kerntype = 'py1vN'

    def prepare_arrays(self):
        Approach.prepare_arrays(self)
        nleads, ndm1 = self.si.nleads, self.si.ndm1
        self.phi1fct = np.zeros((nleads, ndm1, 2), dtype=complexnp)
        self.phi1fct_energy = np.zeros((nleads, ndm1, 2), dtype=complexnp)
        self.phi1 = np.zeros((nleads, ndm1), dtype=complexnp)

    def clean_arrays(self):
        Approach.clean_arrays(self)
        self.phi1fct.fill(0.0)
        self.phi1fct_energy.fill(0.0)
        self.phi1.fill(0.0)

    def generate_fct(self):
        """
        Make factors used for generating 1vN, Redfield master equation kernels.

        Parameters
        ----------
        phi1fct : array
            (Modifies) Factors used for generating 1vN, Redfield master equation kernels.
        phi1fct_energy : array
            (Modifies) Factors used to calculate energy and heat currents in 1vN, Redfield approaches.
        """
        E, si = self.qd.Ea, self.si,
        mulst, tlst, dlst = self.leads.mulst, self.leads.tlst, self.leads.dlst
        itype, limit = self.funcp.itype, self.funcp.dqawc_limit
        ncharge, nleads, statesdm = si.ncharge, si.nleads, si.statesdm

        phi1fct = self.phi1fct
        phi1fct_energy = self.phi1fct_energy
        for charge in range(ncharge-1):
            ccharge = charge+1
            bcharge = charge
            for c, b in itertools.product(statesdm[ccharge], statesdm[bcharge]):
                cb = si.get_ind_dm1(c, b, bcharge)
                Ecb = E[c]-E[b]
                for l in range(nleads):
                    rez = func_1vN(Ecb, mulst[l], tlst[l], dlst[l, 0], dlst[l, 1], itype, limit)
                    phi1fct[l, cb, 0] = rez[0]
                    phi1fct[l, cb, 1] = rez[1]
                    phi1fct_energy[l, cb, 0] = rez[2]
                    phi1fct_energy[l, cb, 1] = rez[3]

    def generate_coupling_terms(self, b, bp, bcharge):
        Tba, phi1fct = self.leads.Tba, self.phi1fct
        si, kh = self.si, self.kernel_handler
        nleads, statesdm = si.nleads, si.statesdm

        acharge = bcharge-1
        ccharge = bcharge+1

        # --------------------------------------------------
        for a, ap in itertools.product(statesdm[acharge], statesdm[acharge]):
            if kh.is_included(a, ap, acharge):
                bpa = si.get_ind_dm1(bp, a, acharge)
                bap = si.get_ind_dm1(b, ap, acharge)
                fct_aap = 0
                for l in range(nleads):
                    fct_aap += (+ Tba[l, b, a]*Tba[l, ap, bp]*phi1fct[l, bpa, 0].conjugate()
                                - Tba[l, b, a]*Tba[l, ap, bp]*phi1fct[l, bap, 0])
                kh.set_matrix_element(fct_aap, b, bp, bcharge, a, ap, acharge)
        # --------------------------------------------------
        for bpp in statesdm[bcharge]:
            if kh.is_included(bpp, bp, bcharge):
                fct_bppbp = 0
                for a in statesdm[acharge]:
                    bpa = si.get_ind_dm1(bp, a, acharge)
                    for l in range(nleads):
                        fct_bppbp += +Tba[l, b, a]*Tba[l, a, bpp]*phi1fct[l, bpa, 1].conjugate()
                for c in statesdm[ccharge]:
                    cbp = si.get_ind_dm1(c, bp, bcharge)
                    for l in range(nleads):
                        fct_bppbp += +Tba[l, b, c]*Tba[l, c, bpp]*phi1fct[l, cbp, 0]
                kh.set_matrix_element(fct_bppbp, b, bp, bcharge, bpp, bp, bcharge)
            # --------------------------------------------------
            if kh.is_included(b, bpp, bcharge):
                fct_bbpp = 0
                for a in statesdm[acharge]:
                    ba = si.get_ind_dm1(b, a, acharge)
                    for l in range(nleads):
                        fct_bbpp += -Tba[l, bpp, a]*Tba[l, a, bp]*phi1fct[l, ba, 1]
                for c in statesdm[ccharge]:
                    cb = si.get_ind_dm1(c, b, bcharge)
                    for l in range(nleads):
                        fct_bbpp += -Tba[l, bpp, c]*Tba[l, c, bp]*phi1fct[l, cb, 0].conjugate()
                kh.set_matrix_element(fct_bbpp, b, bp, bcharge, b, bpp, bcharge)
        # --------------------------------------------------
        for c, cp in itertools.product(statesdm[ccharge], statesdm[ccharge]):
            if kh.is_included(c, cp, ccharge):
                cbp = si.get_ind_dm1(c, bp, bcharge)
                cpb = si.get_ind_dm1(cp, b, bcharge)
                fct_ccp = 0
                for l in range(nleads):
                    fct_ccp += (+ Tba[l, b, c]*Tba[l, cp, bp]*phi1fct[l, cbp, 1]
                                - Tba[l, b, c]*Tba[l, cp, bp]*phi1fct[l, cpb, 1].conjugate())
                kh.set_matrix_element(fct_ccp, b, bp, bcharge, c, cp, ccharge)
        # --------------------------------------------------

    def generate_current(self):
        """
        Calculates currents using 1vN approach.

        Parameters
        ----------
        phi1 : array
            (Modifies) Values of first order density matrix elements
            stored in nleads by ndm1 numpy array.
        current : array
            (Modifies) Values of the current having nleads entries.
        energy_current : array
            (Modifies) Values of the energy current having nleads entries.
        heat_current : array
            (Modifies) Values of the heat current having nleads entries.
        """
        E, Tba = self.qd.Ea, self.leads.Tba
        phi1fct, phi1fct_energy = self.phi1fct, self.phi1fct_energy

        si = self.si
        ncharge, nleads, statesdm = si.ncharge, si.nleads, si.statesdm

        phi1 = self.phi1
        current = self.current
        energy_current = self.energy_current

        kh = self.kernel_handler
        for charge in range(ncharge-1):
            ccharge = charge+1
            bcharge = charge
            for c, b in itertools.product(statesdm[ccharge], statesdm[bcharge]):
                cb = si.get_ind_dm1(c, b, bcharge)

                for l in range(nleads):
                    current_l, energy_current_l, phi1_l = 0, 0, 0

                    fct1 = phi1fct[l, cb, 0]
                    fct2 = phi1fct[l, cb, 1]
                    fct1h = phi1fct_energy[l, cb, 0]
                    fct2h = phi1fct_energy[l, cb, 1]

                    for bp in si.statesdm[bcharge]:
                        if not kh.is_included(bp, b, bcharge):
                            continue
                        phi0bpb = kh.get_phi0_element(bp, b, bcharge)

                        phi1_l += Tba[l, c, bp]*phi0bpb*fct1
                        current_l += Tba[l, b, c]*Tba[l, c, bp]*phi0bpb*fct1
                        energy_current_l += Tba[l, b, c]*Tba[l, c, bp]*phi0bpb*fct1h

                    for cp in statesdm[ccharge]:
                        if not kh.is_included(c, cp, ccharge):
                            continue
                        phi0ccp = kh.get_phi0_element(c, cp, ccharge)

                        phi1_l += Tba[l, cp, b]*phi0ccp*fct2
                        current_l += Tba[l, b, c]*phi0ccp*Tba[l, cp, b]*fct2
                        energy_current_l += Tba[l, b, c]*phi0ccp*Tba[l, cp, b]*fct2h

                    phi1[l, cb] += phi1_l
                    current[l] += -2*current_l.imag
                    energy_current[l] += -2*energy_current_l.imag

        self.heat_current[:] = energy_current - current*self.leads.mulst class Approach2vN(ApproachBase2vN):

    kerntype = 'py2vN'

    #region Preparation

    def prepare_kern(self):
        if self.is_zeroth_iteration:
            self.prepare_energy_grid()
            self.prepare_distribution()

        ApproachBase2vN.prepare_kern(self)

    def prepare_energy_grid(self):
        """"Prepares the energy grid."""
        self.make_Ek_grid()

        # Define the extended grid Ek_grid_ext for calculations outside the bandwidth
        # Here funcp.emin, funcp.emax are defined
        self.determine_emin_emax()

        # Here Ek_grid_ext, funcp.kpnt_left, funcp.kpnt_right are defined
        self.determine_grid_ext()

    def make_Ek_grid(self):
        """Make an energy grid on which 2vN equations are solved. """
        if self.funcp.kpnt is None:
            raise ValueError('kpnt needs to be specified.')
        if self.si.nleads > 0:
            dmin = np.min(self.leads.dlst)
            dmax = np.max(self.leads.dlst)
            Ek_grid, kpnt = self.Ek_grid, self.funcp.kpnt
            if Ek_grid[0] != dmin or Ek_grid[-1] != dmax or Ek_grid.shape[0] != kpnt:
                self.funcp.dmin = dmin
                self.funcp.dmax = dmax
                self.Ek_grid = np.linspace(dmin, dmax, kpnt)
                #
                if self.niter != -1:
                    print("WARNING: Ek_grid has changed. Restarting the calculation.")
                    self.restart()
                #
                if ((dmin * np.ones(self.si.nleads)).tolist() != self.leads.dlst.T[0].tolist() or
                    (dmax * np.ones(self.si.nleads)).tolist() != self.leads.dlst.T[1].tolist()):
                    print("WARNING: The bandwidth and Ek_grid for all leads will be the same: from " +
                          "dmin=" + str(dmin) + " to dmax=" + str(dmax) + ".")

    def determine_emin_emax(self):
        """
        Finds how much the Ek_grid has to be expanded
        above the bandwidth D of the leads.

        Parameters
        ----------
        funcp.emin : float
            (Modifies) Minimal energy in the updated Ek_grid.
        funcp.emax : float
            (Modifies) Maximal energy in the updated Ek_grid.
        """
        E, si = self.qd.Ea, self.si
        dmin, dmax = self.funcp.dmin, self.funcp.dmax

        lst = [dmin, dmax]
        for charge in range(si.ncharge):
            for b, bp in itertools.product(si.statesdm[charge], si.statesdm[charge]):
                lst.append(dmax-E[b]+E[bp])
                lst.append(dmin-E[b]+E[bp])
        for charge in range(si.ncharge-2):
            for d, b in itertools.product(si.statesdm[charge+2], si.statesdm[charge]):
                lst.append(dmin+E[d]-E[b])
                lst.append(dmax+E[d]-E[b])

        self.funcp.emax = max(lst)
        self.funcp.emin = min(lst)

    def determine_grid_ext(self):
        """
        Expands the Ek_grid from emin to emax.

        Parameters
        ----------
        Ek_grid_ext : ndarray
            (Modifies) Extended Ek_grid from emin to emax.
        funcp.kpnt_left : int
            (Modifies) Number of points Ek_grid is extended to the left.
        funcp.kpnt_right : int
            (Modifies) Number of points Ek_grid is extended to the right.
        """
        emin_, emax_ = self.funcp.emin, self.funcp.emax
        dmin, dmax = self.funcp.dmin, self.funcp.dmax
        ext_fct = self.funcp.ext_fct
        Ek_grid = self.Ek_grid

        step = Ek_grid[1]-Ek_grid[0]
        emin = ext_fct*(emin_-dmin)+dmin
        emax = ext_fct*(emax_-dmax)+dmax
        ext_left = np.sort(-np.arange(-dmin+step, -emin+step, step))
        ext_right = np.arange(dmax+step, emax+step, step)

        self.Ek_grid_ext = np.concatenate((ext_left, Ek_grid, ext_right))
        self.funcp.kpnt_left, self.funcp.kpnt_right = len(ext_left), len(ext_right)

    def prepare_distribution(self):
        """"
        Prepares the Fermi distribution functions
        and their Hilbert transform on the energy grid.

        Parameters
        ----------
        fkp : ndarray
            (Modifies) nleads by len(Ek_grid_ext) numpy array containing Fermi function.
        fkm : ndarray
            (Modifies) nleads by len(Ek_grid_ext) numpy array containing 1-Fermi function.
        hfkp : ndarray
            (Modifies) Hilbert transform of fkp.
        hfkm : ndarray
            (Modifies) Hilbert transform of fkm.
        """
        # Generate the Fermi functions on the grid
        # This is necessary to generate only if Ek_grid, mulst, or tlst are changed
        nleads = self.si.nleads
        mulst, tlst = self.leads.mulst, self.leads.tlst
        Ek_grid = self.Ek_grid
        Eklen = len(Ek_grid)

        self.fkp = np.zeros((nleads, Eklen), dtype=doublenp)
        for l in range(nleads):
            self.fkp[l] = 1/(np.exp((Ek_grid - mulst[l])/tlst[l]) + 1)
        self.fkm = 1-self.fkp
        self.fkp, self.hfkp = get_htransf_fk(self.fkp, self.funcp)
        self.fkm, self.hfkm = get_htransf_fk(self.fkm, self.funcp)

    #endregion Preparation

    #region Generation

    def iterate(self):
        terms_calculator = TermsCalculator2vN(self)
        terms_calculator.iterate()

    def determine_phi1_phi0(self):
        """
        Integrates phi1k over energy.

        Parameters
        ----------
        phi1_phi0 : ndarray
            (Modifies) Numpy array with dimensions (nleads, ndm1, ndm0),
            corresponding to energy integrated Phi[1](k).
        e_phi1_phi0 : ndarray
            (Modifies) Numpy array with dimensions (nleads, ndm1, ndm0),
            corresponding to energy integrated Ek*Phi[1](k).
        """
        phi1k, Ek_grid = self.phi1k, self.Ek_grid
        # Get integrated Phi[1]_{cb} in terms of Phi[0]_{bb'}
        phi1_phi0 = self.phi1_phi0
        e_phi1_phi0 = self.e_phi1_phi0

        Eklen = len(Ek_grid)
        for i in range(Eklen):
            # Trapezoidal rule
            if i == 0:
                dx = Ek_grid[i+1] - Ek_grid[i]
            elif i == Eklen-1:
                dx = Ek_grid[i] - Ek_grid[i-1]
            else:
                dx = Ek_grid[i+1] - Ek_grid[i-1]
            phi1_phi0 += 1/2*dx*phi1k[i]
            e_phi1_phi0 += 1/2*dx*Ek_grid[i]*phi1k[i]

    def generate_kern(self):
        """
        From Phi[1](k) generate equations containing just Phi[0].

        Parameters
        ----------
        kern : ndarray
            (Modifies) Kernel for master equation involving just Phi[0].
        bvec : ndarray
            (Modifies) Right hand side column vector for master equation.
        """
        phi1_phi0, E, Tba, si = self.phi1_phi0, self.qd.Ea, self.leads.Tba, self.si
        ncharge, nleads, ndm0, statesdm = si.ncharge, si.nleads, si.ndm0, si.statesdm

        # Integrated Phi[1]_{bc} in terms of phi1_phi0
        shuffle = np.zeros((ndm0, ndm0), dtype=intnp)
        for charge in range(ncharge):
            for b, bp in itertools.combinations_with_replacement(statesdm[charge], 2):
                bbp = si.get_ind_dm0(b, bp, charge)
                bpb = si.get_ind_dm0(bp, b, charge)
                shuffle[bbp, bpb] = 1
                shuffle[bpb, bbp] = 1
        phi1_phi0_conj = np.dot(np.conjugate(phi1_phi0), shuffle)

        kern = self.kern
        for charge in range(ncharge):
            acharge = charge-1
            bcharge = charge
            ccharge = charge+1
            for b, bp in itertools.product(statesdm[bcharge], statesdm[bcharge]):
                bbp = si.get_ind_dm0(b, bp, bcharge)
                kern[bbp, bbp] = E[b]-E[bp]

                for a1 in statesdm[acharge]:
                    bpa1 = si.get_ind_dm1(bp, a1, acharge)
                    ba1 = si.get_ind_dm1(b, a1, acharge)
                    for l in range(nleads):
                        kern[bbp] += +Tba[l, b, a1]*phi1_phi0_conj[l, bpa1]
                        kern[bbp] += -phi1_phi0[l, ba1]*Tba[l, a1, bp]

                for c1 in statesdm[ccharge]:
                    c1bp = si.get_ind_dm1(c1, bp, bcharge)
                    c1b = si.get_ind_dm1(c1, b, bcharge)
                    for l in range(nleads):
                        kern[bbp] += +Tba[l, b, c1]*phi1_phi0[l, c1bp]
                        kern[bbp] += -phi1_phi0_conj[l, c1b]*Tba[l, c1, bp]

    def generate_current(self):
        """
        Finds the currents for 2vN approach after an approximate
        solution for integral 2vN equations was found.

        Parameters
        ----------
        phi1 : ndarray
            (Modifies) nleads by ndm1 numpy array, which gives energy integrated
            Phi[1](k) after the values for Phi[0] were inserted.
        current : ndarray
            (Modifies) Values of the current having nleads entries.
        energy_current : ndarray
            (Modifies) Values of the energy current having nleads entries.
        heat_current : ndarray
            (Modifies) Values of the heat current having nleads entries.
        """
        phi1_phi0, e_phi1_phi0 = self.phi1_phi0, self.e_phi1_phi0
        phi0, Tba, si = self.phi0, self.leads.Tba, self.si
        ncharge, nleads, ndm1, statesdm = si.ncharge, si.nleads, si.ndm1, si.statesdm

        phi1 = self.phi1
        h_phi1 = self.h_phi1
        current = self.current
        energy_current = self.energy_current

        for l in range(nleads):
            phi1[l] = np.dot(phi1_phi0[l], phi0)
            h_phi1[l] = np.dot(e_phi1_phi0[l], phi0)

        for charge in range(ncharge-1):
            ccharge = charge+1
            bcharge = charge
            for c, b in itertools.product(statesdm[ccharge], statesdm[bcharge]):
                cb = si.get_ind_dm1(c, b, bcharge)

                for l in range(nleads):
                    current_l = phi1[l, cb]*Tba[l, b, c]
                    energy_current_l = h_phi1[l, cb]*Tba[l, b, c]

                    current[l] += -2*current_l.imag
                    energy_current[l] += -2*energy_current_l.imag

        self.heat_current[:] = energy_current - current*self.leads.mulst

    #endregion Generation class ApproachPyRTD(Approach):

    kerntype = 'pyRTD'

    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)
        self.BW_Ozaki_expansion = 0
        self.Ozaki_poles_and_residues = None
        self.off_diag_corrections = self.funcp.off_diag_corrections
        self.ImGamma = False
        self.printed_warning_ImGamma = False
        self.nsingle_warning_printed = False

    def get_kern_size(self):
        return self.si.npauli

    def restart(self):
        Approach.restart(self)
        self.Wdd = None
        self.WE1 = None
        self.WE2 = None
        self.ReWnd = None
        self.ImWnd = None
        self.ReWdn = None
        self.ImWdn = None
        self.Lnn = None

    def prepare_kernel_handler(self):
        self.kernel_handler = KernelHandlerRTD(self.si)

    def prepare_arrays(self):
        Approach.prepare_arrays(self)

        nleads, ndm1 = self.si.nleads, self.si.ndm1
        self.paulifct = np.zeros((nleads, ndm1, 2), dtype=doublenp)

        kern_size = self.get_kern_size()
        nleads = self.si.nleads
        self.Wdd = np.zeros((nleads, kern_size, kern_size), dtype=self.dtype)
        self.WE1 = np.zeros((nleads, kern_size, kern_size), dtype=self.dtype)
        self.WE2 = np.zeros((nleads, kern_size, kern_size), dtype=self.dtype)

        kh = self.kernel_handler
        kh.Wdd = self.Wdd
        kh.WE1 = self.WE1
        kh.WE2 = self.WE2

        self.W2 = np.zeros((nleads, kern_size, kern_size), dtype=self.dtype)

        self.generate_LN()
        self.LE = np.zeros((kern_size, kern_size), dtype=self.dtype)

        if self.off_diag_corrections:
            kern_size2 = 2 * self.si.ndm0 - 2 * self.si.npauli
            self.ReWnd = np.zeros((nleads, kern_size2, kern_size), dtype=self.dtype)
            self.ImWnd = np.zeros((nleads, kern_size2, kern_size), dtype=self.dtype)
            self.ReWdn = np.zeros((nleads, kern_size, kern_size2), dtype=self.dtype)
            self.ImWdn = np.zeros((nleads, kern_size, kern_size2), dtype=self.dtype)
            self.Lnn = np.zeros((kern_size2, kern_size2), dtype=self.dtype)

            kh.ReWnd = self.ReWnd
            kh.ImWnd = self.ImWnd
            kh.ReWdn = self.ReWdn
            kh.ImWdn = self.ImWdn
            kh.Lnn = self.Lnn

        kh.set_matrix_list()

    def clean_arrays(self):
        Approach.clean_arrays(self)
        self.ImGamma = False
        self.paulifct.fill(0.0)

        self.Wdd.fill(0.0)
        self.WE1.fill(0.0)
        self.WE2.fill(0.0)
        self.W2.fill(0.0)
        self.LE.fill(0.0)

        if self.off_diag_corrections:
            self.ReWnd.fill(0.0)
            self.ImWnd.fill(0.0)
            self.ReWdn.fill(0.0)
            self.ImWdn.fill(0.0)
            self.Lnn.fill(0.0)

    def generate_kern(self):
        r""" Generates all kernels including tunnel processes of orders :math:`t^2` and :math:`t^4`.

        The total kernel used to solve for :math:`\phi_0` is :math:`W =\sum_r W^r= W_{dd}^{(1)} + W_{dd}^{(2)}
        + W_{dn}^{(1)} (L_{nn})^{-1} W_{nd}^{(1)}`. The last term is ignored if `off_diag_corrections` is False.

        Parameters
        ----------
        self.kern : ndarray
            (Modifies) The total Kernel for the diagonal density matrix. Has npauli * npauli entries.
        self.Wdd :  ndarray
            (Modifies) The lead-resolved Kernel for the diagonal density matrix. Has
            nleads * npauli * npauli entries.

        """
        si, kh = self.si, self.kernel_handler
        ncharge, statesdm = si.ncharge, si.statesdm
        self.off_diag_corrections = self.funcp.off_diag_corrections

        if (not np.all(np.isclose(self.leads.tlst, self.leads.tlst[0]))) or np.any(abs(self.leads.Tba.imag)>0):
            self.set_Ozaki_params()

        for bcharge in range(ncharge):
            for b in statesdm[bcharge]:
                if not kh.is_unique(b, b, bcharge):
                    continue
                self.generate_row_1st_order_kernel(b, bcharge)
                self.generate_col_diag_kern_2nd_order(b, bcharge)
                self.generate_row_1st_energy_kernel(b, bcharge)
                self.generate_row_2nd_energy_kernel(b, bcharge)

                if self.off_diag_corrections:
                    self.generate_col_nondiag_kern_1st_order_nd(b, bcharge)

        kern_size = self.get_kern_size()
        self.kern[:kern_size, :kern_size] += np.sum(self.Wdd, 0)

        if self.off_diag_corrections:
            for bcharge in range(ncharge):
                for b in statesdm[bcharge]:
                    for bp in statesdm[bcharge]:
                        if b == bp:
                            continue
                        self.generate_col_nondiag_kern_1st_order_dn(b, bp, bcharge)
                        self.generate_row_inverse_Liouvillian(b, bp, bcharge)
            self.add_off_diag_corrections()

    def add_off_diag_corrections(self):
        """
        Adds :math:`W_{dn}^{(1)} (L_{nn})^{-1} W_{nd}^{(1)}` to the total kernel.

        Parameters
        ----------
        self.kern : ndarray
            (Modifies) The total kernel

        self.Wdd : ndarray
            (Modifies) The lead-resolved diagonal kernel
        """
        kh = self.kernel_handler
        Wcorr = np.zeros(self.Wdd.shape, dtype=doublenp)
        for l in range(self.si.nleads):
            # Since off-diagonal kernels contain imaginary numbers we get two contributions to
            # the kernel
            Wcorr[l, :, :] += np.matmul(np.matmul(kh.ReWdn[l, :, :], kh.Lnn[:, :]),
                                        np.sum(kh.ImWnd[:, :], 0))
            Wcorr[l, :, :] += np.matmul(np.matmul(kh.ImWdn[l, :, :], kh.Lnn[:, :]),
                                        np.sum(kh.ReWnd[:, :], 0))
        self.Wdd += Wcorr
        kern_size = self.get_kern_size()
        self.kern[:kern_size, :kern_size] += np.sum(Wcorr, 0)

    def generate_LN(self):
        """Generates the diagonal of :math:`L_{N_{dot}}^+`

        Parameters
        ----------
        self.LN : ndarray
            (Modifies) the digononal entries of the anti-commutator Liouvillian for the number
            operator. The array has npauli entries.
        """
        charge_lst = []
        for charge in range(self.si.ncharge):
            for a in self.si.statesdm[charge]:
                if self.si.get_ind_dm0(a, a, charge, 2):
                    charge_lst.append(2 * charge)
        self.LN = np.array(charge_lst)

    def generate_current(self):
        r""" Calculates currents for the RTD approach.

        Charge current for reservoir r is evaluated as :math:`I^r = 1/2 \cdot Tr (L_{N_{dot}}^+ W^r \phi_0)`

        Energy current for reservoir r is evaluated as :math:`E^r = 1/2\cdot Tr (L_H^+ W^r \phi_0) - 1/2\cdot
        Tr(W_{E,1}^r \phi_0) + 1/2\cdot Tr(W_{E,2}^r \phi_0)`

        Heat current for reservoir r is evaluated as :math:`Q^r = E^r - \mu_r\cdot I^r`.

        Parameters
        ----------
        self.current : array
            (Modifies) Charge current in each lead.
        self.energy_current : array
            (Modifies) Energy charge current in each lead.
        self.heat_current : array
            (Modifies) Heat charge current in each lead.

        """
        kh = self.kernel_handler
        nleads = self.si.nleads
        LE = 2 * self.qd.Ea

        current = np.zeros(nleads, dtype=doublenp)
        energy_current = np.zeros(nleads, dtype=doublenp)

        for l in range(nleads):
            current[l] = 0.5 * np.dot(self.LN, np.dot(kh.Wdd[l, :, :], self.phi0))
            energy_current[l] = 0.5 * np.dot(LE, np.dot(kh.Wdd[l, :, :], self.phi0))
            energy_current[l] += -0.5 * np.sum(np.dot(kh.WE1[l, :, :], self.phi0))
            energy_current[l] += 0.5 * np.sum(np.dot(kh.WE2[l, :, :], self.phi0))

        self.current = current
        self.energy_current = energy_current
        self.heat_current = energy_current - current * self.leads.mulst

        if self.ImGamma:
            self.energy_current.fill(np.nan)
            self.heat_current.fill(np.nan)
            if not self.printed_warning_ImGamma:
                print('Warning! Complex matrix elements detected, which are not supported for the RTD approach ' +
                      'when calculating the energy current.')
                self.printed_warning_ImGamma = True

        if self.si.nsingle == 0:
            if not self.nsingle_warning_printed:
                print('Warning! No single particle tunneling amplitudes (tleads) detected. Corrections to the energy ' +
                      'current in the RTD approach uses tleads. Please specify BuilderManyBody.tleads_array and ' +
                      'BuilderManyBody.nsingle, if possible.\n\nThe correction terms can be neglected if no single' +
                      ' particle state is connected to more than one lead.')
                self.nsingle_warning_printed = True

    def generate_fct(self):
        """
        Make factors used for generating the first order diagonal kernel :math:`W_{dd}^{(1)}`.

        Parameters
        ----------
        paulifct : array
            (Modifies) Factors used for generating Pauli master equation kernel.
        """
        E, Tba, si = self.qd.Ea, self.leads.Tba, self.si
        mulst, tlst, dlst = self.leads.mulst, self.leads.tlst, self.leads.dlst
        ncharge, nleads, statesdm = si.ncharge, si.nleads, si.statesdm

        itype = self.funcp.itype
        paulifct = self.paulifct
        for charge in range(ncharge-1):
            ccharge = charge+1
            bcharge = charge
            for c, b in itertools.product(statesdm[ccharge], statesdm[bcharge]):
                cb = si.get_ind_dm1(c, b, bcharge)
                Ecb = E[c]-E[b]
                for l in range(nleads):
                    xcb = (Tba[l, b, c]*Tba[l, c, b]).real
                    rez = func_pauli(Ecb, mulst[l], tlst[l], dlst[l, 0], dlst[l, 1], itype)
                    paulifct[l, cb, 0] = xcb*rez[0]
                    paulifct[l, cb, 1] = xcb*rez[1]

    def generate_row_1st_order_kernel(self, b, bcharge):
        """Generates a row in the first order diagonal kernel :math:`W_{dd}^{(1)}`.

        Parameters
        ----------
        b : int
            the state (row)

        bcharge : int
            charge of state b

        self.Wdd : ndarray
            (Modifies) The kernel connecting diagional density-matrix elements. This Kernel
            has npauli * npauli entries.
        """
        paulifct = self.paulifct
        si, kh = self.si, self.kernel_handler
        nleads, statesdm = si.nleads, si.statesdm

        acharge = bcharge-1
        ccharge = bcharge+1

        bb = si.get_ind_dm0(b, b, bcharge)
        for a in statesdm[acharge]:
            aa = si.get_ind_dm0(a, a, acharge)
            ba = si.get_ind_dm1(b, a, acharge)
            for l in range(nleads):
                fctm = -paulifct[l, ba, 1]
                fctp = paulifct[l, ba, 0]
                kh.set_matrix_element_dd(l, fctm, fctp, bb, aa, 0)
        for c in statesdm[ccharge]:
            cc = si.get_ind_dm0(c, c, ccharge)
            cb = si.get_ind_dm1(c, b, bcharge)
            for l in range(nleads):
                fctm = -paulifct[l, cb, 0]
                fctp = paulifct[l, cb, 1]
                kh.set_matrix_element_dd(l, fctm, fctp, bb, cc, 0)

    def generate_row_1st_energy_kernel(self, b, bcharge):
        r""" Generates a row of the first kernel for the barrier part of the energy current :math:`W_{E,1}`. This kernel
         is obtained from a diagram of the form
         :math:`L_{T,r}^+  L_{T, r_2}^-  (z-L_{dot} - L_r )^{-1}  L_T  (z-L_{dot} - L_r )^{-1} L_T \phi_0`,
         when contracting the first and third, as well as second and fourth, reservoir superoperators. Assumes
         that the wide-band limit is valid and that all products of tunnel matrix elements are real.

         Parameters
         ----------
         b : int
            the row (state)
         charge : int
            the charge of state b
         self.W2E_1 : ndarray
             (Modifies) The second order kernel used for evaluating the energy current. This Kernel
             has npauli * npauli entries.
         """
        (E, Tba, tleads, mulst, tlst, dlst) = (self.qd.Ea, self.leads.Tba, self.leads.tleads_array,
                                                   self.leads.mulst, self.leads.tlst, self.leads.dlst)
        si, kh = self.si, self.kernel_handler
        nleads, statesdm, nsingle = si.nleads, si.statesdm, si.nsingle

        acharge = bcharge-1
        ccharge = bcharge+1
        maxTemp = max(tlst)
        t_cutoff = 1e-15 * maxTemp ** 2

        bb = si.get_ind_dm0(b, b, bcharge)
        for a in statesdm[acharge]:
            aa = si.get_ind_dm0(a, a, acharge)
            for l in range(nleads):
                mu, Tr, gamma = mulst[l], tlst[l], 0.0
                dE = E[b] - E[a]
                for lp in range(nleads):
                    if lp == l: continue
                    for n1 in range(nsingle):
                        gamma += Tba[l, a, b] * Tba[lp, a, b].conj() * tleads[l, n1] * tleads[lp, n1].conj()
                temp = gamma.real * phi((dE - mu) / Tr, dlst[l, 0] / Tr, dlst[l, 1] / Tr)
                temp += gamma.real * phi(-(dE - mu) / Tr, dlst[l, 0] / Tr, dlst[l, 1] / Tr)
                if abs(gamma.imag) > t_cutoff:
                    self.ImGamma = True
                    #temp += gamma.imag * fermi_func((dE - mu) / Tr)*np.pi
                    #temp += gamma.imag * fermi_func(-(dE - mu) / Tr) * np.pi

                temp *= np.pi
                kh.set_matrix_element_dd(l, temp, temp, bb, aa, 1)

        for c in statesdm[ccharge]:
            cc = si.get_ind_dm0(c, c, ccharge)
            for l in range(nleads):
                mu, Tr, gamma = mulst[l], tlst[l], 0.0
                dE = E[c] - E[b]
                for lp in range(nleads):
                    if lp == l: continue
                    for n1 in range(nsingle):
                        gamma += Tba[l, b, c] * Tba[lp, b, c].conj() * tleads[l, n1] * tleads[lp, n1].conj()
                temp = gamma.real * phi((dE - mu) / Tr, dlst[l, 0] / Tr, dlst[l, 1] / Tr)
                temp += gamma.real * phi(-(dE - mu) / Tr, dlst[l, 0] / Tr, dlst[l, 1] / Tr)
                if abs(gamma.imag) > t_cutoff:
                    self.ImGamma = True
                    #temp += gamma.imag * fermi_func((dE - mu) / Tr) * np.pi
                    #temp += gamma.imag * fermi_func(-(dE - mu) / Tr) * np.pi

                temp *= np.pi
                kh.set_matrix_element_dd(l, temp, temp, bb, cc, 1)

    def generate_row_2nd_energy_kernel(self, b, bcharge):
        r""" Generates a row in the second kernel for the barrier part of the energy current :math:`W_{E,2}`.
        This kernel is obatined from a diagram of the form
        :math:`L_{T,r}^+  L_{T, r_2}^-  (z-L_{dot} - L_r )^{-1}  L_T  (z-L_{dot} - L_r )^{-1} L_T \phi_0`,
        when contracting the first and fourth, as well as second and third, reservoir superoperators. Assumes
        that the wide-band limit is valid and that all products of tunnel matrix elements are real.

        Parameters
        ----------
        b : int
            the row (state)
        charge : int
            the charge of state b
        self.W2E_2 : ndarray
            (Modifies) The second order kernel used for evaluating the energy current. This Kernel
            has npauli * npauli entries.
        """
        (E, si, Tba, tleads, mulst, tlst, dlst) = (self.qd.Ea, self.si, self.leads.Tba, self.leads.tleads_array,
                                                   self.leads.mulst, self.leads.tlst, self.leads.dlst)
        si, kh = self.si, self.kernel_handler
        nleads, statesdm, nsingle = si.nleads, si.statesdm, si.nsingle

        acharge = bcharge - 1
        ccharge = bcharge + 1
        maxTemp = max(tlst)
        t_cutoff = 1e-15*maxTemp**2

        bb = si.get_ind_dm0(b, b, bcharge)
        for a in statesdm[acharge]:
            aa = si.get_ind_dm0(a, a, acharge)
            for l in range(nleads):
                temp = 0.0
                for lp in range(nleads):
                    if lp != l:
                        mu, Tr, gamma = mulst[lp], tlst[lp], 0.0
                        for n1 in range(nsingle):
                            gamma += Tba[l, a, b] * Tba[lp, a, b].conj() * tleads[l, n1] * tleads[lp, n1].conj()
                        if abs(gamma.imag) > t_cutoff:
                            self.ImGamma = True
                        dE = E[b] - E[a]
                        temp += gamma.real * phi((dE - mu) / Tr, dlst[lp, 0] / Tr, dlst[lp, 1] / Tr)
                        temp += gamma.real * phi(-(dE - mu) / Tr, dlst[lp, 0] / Tr, dlst[lp, 1] / Tr)
                        if abs(gamma.imag) > t_cutoff:
                            self.ImGamma = True
                            #temp += gamma.imag * fermi_func((dE - mu) / Tr) * np.pi
                            #temp += gamma.imag * fermi_func(-(dE - mu) / Tr) * np.pi
                temp *= np.pi
                kh.set_matrix_element_dd(l, temp, temp, bb, aa, 2)

        for c in statesdm[ccharge]:
            cc = si.get_ind_dm0(c, c, ccharge)
            for l in range(nleads):
                temp = 0.0
                for lp in range(nleads):
                    if lp != l:
                        mu, Tr, gamma = mulst[lp], tlst[lp], 0.0
                        for n1 in range(nsingle):
                            gamma += Tba[l, b, c] * Tba[lp, b, c].conj() * tleads[l, n1] * tleads[lp, n1].conj()
                        dE = E[c] - E[b]
                        temp += gamma.real * phi((dE - mu) / Tr, dlst[lp, 0] / Tr, dlst[lp, 1] / Tr)
                        temp += gamma.real * phi(-(dE - mu) / Tr, dlst[lp, 0] / Tr, dlst[lp, 1] / Tr)
                        if abs(gamma.imag) > t_cutoff:
                            self.ImGamma = True
                            #temp += gamma.imag * fermi_func((dE - mu) / Tr) * np.pi
                            #temp += gamma.imag * fermi_func(-(dE - mu) / Tr) * np.pi
                temp *= np.pi
                kh.set_matrix_element_dd(l, temp, temp, bb, cc, 2)

    def generate_col_diag_kern_2nd_order(self, a0, charge):
        """Partly generates a column in the second order kernel for the diagonal density matrix :math:`W_{dd}^{(2)}`.
        Due to symmetries among the diagrammatic contributions for different matrix elements also contributions to
        other columns are generated. Assumes that the wide band limit is valid.

        Parameters
        ----------
        a0 : int
            initial state. Sets the column

        charge : int
            charge of state a0

        self.Wdd : ndarray
            (Modifies) diagonal lead-resolved kernel.

        """
        # Evaluating the matrix elements requires summing over three pairs of states |a_+><a_-|,
        # two leads, two electron-hole indices and four propagator signs. This is done as follows:
        # one starts with a loop over the initial states. Then, nested loops iterate over the
        # possible intermediate states of the diagram (which is a fairly small number due
        # to several delta functions arising when evaluating the diagrams). Eventually, the final state
        # is reached and a contribution to the matrix element, specified by the initial and the final
        # state, is found (this only provides a contribution since the intermediate states were not
        # fully looped over yet.) Calculating the full matrix elements requires calling this function
        # for every column in the kernel.
        #
        # Symmetries between diagram contributions are used to avoid looping over some Keldysh- and
        # electron-hole indices. Specifically p1 = 1, p4 = 1 and eta1 = 1 are fixed.
        #
        # For the tunnel matrix elements the following rules apply:
        # 1) t = Tba(i->f) if charge_i < charge_f else t = Tba(f->i).conj()
        # 2) t_n = t_n.conj() if p_n == -1
        #
        # Variable names follow Leijnse et al PRB 78, 235424 (2008). Specifically:
        # - aNp/aNm: states
        # - eta : electron-hole index
        # - r : lead index
        # - p : Keldysh sign
        # - z : energy difference

        statesdm, Tba, E = self.si.statesdm, self.leads.Tba, self.qd.Ea
        tlst, mulst, dlst = self.leads.tlst, self.leads.mulst, self.leads.dlst
        kh = self.kernel_handler
        nleads = self.si.nleads
        b_and_R = self.Ozaki_poles_and_residues

        t_cutoff1 = 0.0
        t_cutoff2 = 1e-10*max(tlst)
        t_cutoff3 = 1e-20*max(tlst)**2
        indx0 = self.si.get_ind_dm0(a0, a0, charge)
        for r0, r1 in product(range(nleads), range(nleads)):
            T1, T2 = tlst[r0], tlst[r1]
            mu1, mu2 = mulst[r0], mulst[r1]
            D = np.abs(dlst[r0, 1]) + np.abs(dlst[r0, 0])
            #N1 = (N0, N0 + 1), a1- = a0
            for a1p in statesdm[charge+1]:
                t = Tba[r0, a1p, a0]
                if abs(t) == t_cutoff1:
                    continue
                indx1 = self.si.get_ind_dm0(a1p, a1p, charge + 1)
                E1 = E[a1p] - E[a0]
                #eta1 = 1
                #p1 = 1
                #N2 = (N0, N0 + 2), a2m = a0
                for a2p in statesdm[charge+2]:
                    #p2 = 1
                    t1 = t * Tba[r1, a2p, a1p]
                    if abs(t1) < t_cutoff2:
                        continue
                    E2 = E[a2p] - E[a0]
                    #3 = (N0, N0 + 1 ), a3- = a2-
                    for a3p in statesdm[charge+1]:
                        #charge4 = charge + 0, a4 = a0
                        t2D = t1 * Tba[r1, a2p, a3p].conj() * Tba[r0, a3p, a0].conj()
                        t2X = t1 * Tba[r0, a2p, a3p].conj() * Tba[r1, a3p, a0].conj()
                        E3 = E[a3p] - E[a0]
                        if abs(t2D) > t_cutoff3:
                            tempD = t2D * integralD(1, 1, E1, E2, E3, T1, T2, mu1, mu2, D, b_and_R, abs(t2D.imag)>t_cutoff3)
                            kh.add_element_2nd_order(r0, tempD.real, indx0, indx1, a3p, charge + 1, a0, charge)
                        if abs(t2X) > t_cutoff3:
                            tempX = -t2X * integralX(1, 1, E1, E2, E3, T1, T2, mu1, mu2, D, b_and_R, abs(t2X.imag)>t_cutoff3)
                            kh.add_element_2nd_order(r1, tempX.real, indx0, indx1, a3p, charge + 1, a0, charge)
                    #p2 = -1
                    #N3 = ( N0 +1, N0 + 2)
                    for a3m in statesdm[charge+1]:
                        #charge4 = charge + 1, a4 = a3m
                        t2D = t1 * Tba[r1, a3m, a0].conj() * Tba[r0, a2p, a3m].conj()
                        t2X = t1 * Tba[r0, a3m, a0].conj() * Tba[r1, a2p, a3m].conj()
                        E3 = E[a2p] - E[a3m]
                        if abs(t2D) > t_cutoff3:
                            tempD = t2D * integralD(1, 1, E1, E2, E3, T1, T2, mu1, mu2, D, b_and_R, abs(t2D.imag)>t_cutoff3)
                            kh.add_element_2nd_order(r0, tempD.real, indx0, indx1, a2p, charge + 2, a3m, charge + 1)
                        if abs(t2X) > t_cutoff3:
                            tempX = -t2X * integralX(1, 1, E1, E2, E3, T1, T2, mu1, mu2, D, b_and_R, abs(t2X.imag)>t_cutoff3)
                            kh.add_element_2nd_order(r1, tempX.real, indx0, indx1, a2p, charge + 2, a3m, charge+1)
                #p1 = -1
                #N2 = ( N0 - 1, N0 + 1 ), a2+ = a1+
                for a2m in statesdm[charge-1]:
                    t1 = t * Tba[r1, a0, a2m]
                    if abs(t1) < t_cutoff2:
                        continue
                    E2 = E[a1p] - E[a2m]
                    #p2 = 1
                    #N3 = ( N0 - 1 , N0 ), a3- = a2-
                    for a3p in statesdm[charge]:
                        #charge4 = charge - 1, a0 = a2m
                        t2D = t1 * Tba[r1, a1p, a3p].conj() * Tba[r0, a3p, a2m].conj()
                        t2X = t1 * Tba[r0, a1p, a3p].conj() * Tba[r1, a3p, a2m].conj()
                        E3 = E[a3p] - E[a2m]
                        if abs(t2D) > t_cutoff3:
                            tempD = t2D * integralD(-1, 1, E1, E2, E3, T1, T2, mu1, mu2, D, b_and_R, abs(t2D.imag)>t_cutoff3)
                            kh.add_element_2nd_order(r0, tempD.real, indx0, indx1, a3p, charge, a2m, charge - 1)
                        if abs(t2X) > t_cutoff3:
                            tempX = -t2X * integralX(-1, 1, E1, E2, E3, T1, T2, mu1, mu2, D, b_and_R, abs(t2X.imag)>t_cutoff3)
                            kh.add_element_2nd_order(r1, tempX.real, indx0, indx1, a3p, charge, a2m, charge-1)
                    #p2 = -1
                    #N3 = ( N0 , N0 + 1), a3+ = a2+
                    for a3m in statesdm[charge]:
                        #charge4 = charge + 0, a4 = a3m
                        t2D = t1 * Tba[r1, a3m, a2m].conj() * Tba[r0, a1p, a3m].conj()
                        t2X = t1 * Tba[r0, a3m, a2m].conj() * Tba[r1, a1p, a3m].conj()
                        E3 = E[a1p] - E[a3m]
                        if abs(t2D) > t_cutoff3:
                            tempD = t2D * integralD(-1, 1, E1, E2, E3, T1, T2, mu1, mu2, D, b_and_R, abs(t2D.imag)>t_cutoff3)
                            kh.add_element_2nd_order(r0, tempD.real, indx0, indx1, a1p, charge + 1, a3m, charge)
                        if abs(t2X) > t_cutoff3:
                            tempX = -t2X * integralX(-1, 1, E1, E2, E3, T1, T2, mu1, mu2, D, b_and_R, abs(t2X.imag)>t_cutoff3)
                            kh.add_element_2nd_order(r1, tempX.real, indx0, indx1, a1p, charge + 1, a3m, charge)
                #eta1 = -1
                #p1 = 1
                #N2 = (N0, N0), a2- = a0
                for a2p in statesdm[charge]:
                    E2 = E[a2p] - E[a0]
                    t1 = t * Tba[r1, a1p, a2p].conj()
                    if abs(t1) < t_cutoff2:
                        continue
                    #p2 = 1
                    #N3 = ( N0, N0 +1), a3- = a0
                    for a3p in statesdm[charge+1]:
                        #charge4 = charge, a4 = a0
                        t2D = t1 * Tba[r1, a3p, a2p] * Tba[r0, a3p, a0].conj()
                        E3 = E[a3p] - E[a0]
                        if abs(t2D) > t_cutoff3:
                            tempD = t2D * integralD(1, -1, E1, E2, E3, T1, T2, mu1, mu2, D, b_and_R, abs(t2D.imag)>t_cutoff3)
                            kh.add_element_2nd_order(r0, tempD.real, indx0, indx1, a3p, charge + 1, a0, charge)
                    #N3 = (N0, N0-1), a3- = a0
                    for a3p in statesdm[charge-1]:
                        #charge4 = charge, a4 = a0
                        t2X = t1 * Tba[r0, a2p, a3p].conj() * Tba[r1, a0, a3p]
                        E3 = E[a3p] - E[a0]
                        if abs(t2X) > t_cutoff3:
                            tempX = -t2X * integralX(1, -1, E1, E2, E3, T1, T2, mu1, mu2, D, b_and_R, abs(t2X.imag)>t_cutoff3)
                            kh.add_element_2nd_order(r1, tempX.real, indx0, indx1, a3p, charge - 1, a0, charge)
                    #p2 = -1
                    #N3 = ( N0-1, N0 ), a3+ = a2+
                    for a3m in statesdm[charge-1]:
                        #charge4 = charge - 1, a4 = a3m
                        t2D = t1 * Tba[r1, a0, a3m] * Tba[r0, a2p, a3m].conj()
                        E3 = E[a2p] - E[a3m]
                        if abs(t2D) > t_cutoff3:
                            tempD = t2D * integralD(1, -1, E1, E2, E3, T1, T2, mu1, mu2, D, b_and_R, abs(t2D.imag)>t_cutoff3)
                            kh.add_element_2nd_order(r0, tempD.real, indx0, indx1, a2p, charge, a3m, charge - 1)
                    #N3 = (N0 + 1, N0)
                    for a3m in statesdm[charge+1]:
                        #charge4 = charge + 1, a4 = a3m
                        t2X = t1 * Tba[r0, a3m, a0].conj() * Tba[r1, a3m, a2p]
                        E3 = E[a2p] - E[a3m]
                        if abs(t2X) > t_cutoff3:
                            tempX = -t2X * integralX(1, -1, E1, E2, E3, T1, T2, mu1, mu2, D, b_and_R, abs(t2X.imag)>t_cutoff3)
                            kh.add_element_2nd_order(r1, tempX.real, indx0, indx1, a2p, charge, a3m, charge + 1)
                #p1 = -1
                #N2 = ( N0 + 1  , N0 + 1), a2+ = a1+
                for a2m in statesdm[charge+1]:
                    E2 = E[a1p] - E[a2m]
                    t1 = t * Tba[r1, a2m, a0].conj()
                    #p2 = 1
                    # N3 = (N0 + 1, N0 + 2), a3- = a2-
                    for a3p in statesdm[charge+2]:
                        #charge4 = charge + 1, a4 = a2m
                        t2D = t1 * Tba[r1, a3p, a1p] * Tba[r0, a3p, a2m].conj()
                        E3 = E[a3p] - E[a2m]
                        if abs(t2D) > t_cutoff3:
                            tempD = t2D * integralD(-1, -1, E1, E2, E3, T1, T2, mu1, mu2, D, b_and_R, abs(t2D.imag)>t_cutoff3)
                            kh.add_element_2nd_order(r0, tempD.real, indx0, indx1, a3p, charge + 2, a2m, charge + 1)
                    #N3 = ( N0 + 1, N0 )
                    for a3p in statesdm[charge]:
                        #charge4 = charge + 1, a4 = a2m
                        t2X = t1 * Tba[r0, a1p, a3p].conj() * Tba[r1, a2m, a3p]
                        E3 = E[a3p] - E[a2m]
                        if abs(t2X) > t_cutoff3:
                            tempX = -t2X * integralX(-1, -1, E1, E2, E3, T1, T2, mu1, mu2, D, b_and_R, abs(t2X.imag)>t_cutoff3)
                            kh.add_element_2nd_order(r1, tempX.real, indx0, indx1, a3p, charge, a2m, charge + 1)
                    #p2 = -1
                    #N3 = ( N0, N0+1), a3+ = a2+
                    for a3m in statesdm[charge]:
                        #charge4 = charge, a4 = a3m
                        t2D = t1 * Tba[r1, a2m, a3m] * Tba[r0, a1p, a3m].conj()
                        E3 = E[a1p] - E[a3m]
                        if abs(t2D) > t_cutoff3:
                            tempD = t2D * integralD(-1, -1, E1, E2, E3, T1, T2, mu1, mu2, D, b_and_R, abs(t2D.imag)>t_cutoff3)
                            kh.add_element_2nd_order(r0, tempD.real, indx0, indx1, a1p, charge + 1, a3m, charge)
                    #N3 = ( N0 + 2, N0 + 1 ), a3+ = a2+
                    for a3m in statesdm[charge+2]:
                        #charge4 = charge + 2, a4 = a3m
                        t2X = t1 * Tba[r0, a3m, a2m].conj() * Tba[r1, a3m, a1p]
                        E3 = E[a1p] - E[a3m]
                        if abs(t2X) > t_cutoff3:
                            tempX = -t2X * integralX(-1, -1, E1, E2, E3, T1, T2, mu1, mu2, D, b_and_R, abs(t2X.imag)>t_cutoff3)
                            kh.add_element_2nd_order(r1, tempX.real, indx0, indx1, a1p, charge + 1, a3m, charge + 2)

    def generate_col_nondiag_kern_1st_order_dn(self, a1, b1, charge):
        r""" Calculates a column in :math:`W_{dn}^{(1)}`, the part of the full first order off-diagonal kernel
        connecting diagonal :math:`|a2><a2|` and non-diagonal :math:`|a1><b1|` entries of the density matrix.
        Separates the real and imaginary parts into two matrices. Assumes that the wide band limit is valid.

        Parameters
        ----------
        a1 : int
            first state indexing the column

        b1 : int
            second state indexing the column

        charge: int
            charge of state a1 and b1

        self.ImWdn : ndarray
            (Modifies) :math:`\Im(W_{dn}^{(1)})`, the imaginary part of the first order kernel that connects
            diagonal and non-diagonal elements of the density matrix.
        self.ReWdn : ndarray
            (Modifies) :math:`\Re(W_{dn}^{(1)})`, the real part of the first order kernel that connects
            diagonal and non-diagonal elements of the density matrix.
        """
        # Variables naming follow Leijnse et al.PRB 78, 235424(2008).

        (E, si, Tba, mulst, tlst, dlst) = (self.qd.Ea, self.si, self.leads.Tba,
                                           self.leads.mulst, self.leads.tlst, self.leads.dlst)
        PI = np.pi
        si, kh = self.si, self.kernel_handler
        nleads, statesdm, ncharge = si.nleads, si.statesdm, si.ncharge

        # final state in higher charge state
        if charge != ncharge - 1:
            # Loop over diagonal final state
            for a2 in statesdm[charge + 1]:
                E1 = E[a2] - E[a1]
                E2 = E[a2] - E[b1]
                for l in range(nleads):
                    t2 = Tba[l, a2, a1] * Tba[l, a2, b1].conj()
                    f = fermi_func((E1 - mulst[l]) / tlst[l]) + fermi_func((E2 - mulst[l]) / tlst[l])
                    phi0 = delta_phi((E1 - mulst[l]) / tlst[l], (E2 - mulst[l]) / tlst[l],
                                     dlst[l, 0] / tlst[l], dlst[l, 1] / tlst[l])
                    temp1 = PI * t2.real * f - t2.imag * phi0
                    temp2 = t2.real * phi0 + PI * t2.imag * f
                    kh.add_matrix_element(temp1, l, a2, a2, charge + 1, a1, b1, charge, 3)
                    kh.add_matrix_element(temp2, l, a2, a2, charge + 1, a1, b1, charge, 4)
        # Final state in lower charge state
        if charge != 0:
            # Loop through diagonal final states
            for a2 in statesdm[charge - 1]:
                E1 = E[b1] - E[a2]
                E2 = E[a1] - E[a2]
                for l in range(nleads):
                    t2 = Tba[l, b1, a2] * Tba[l, a1, a2].conj()
                    f = fermi_func(-(E1 - mulst[l]) / tlst[l]) + fermi_func(-(E2 - mulst[l]) / tlst[l])
                    phi0 = delta_phi((E1 - mulst[l]) / tlst[l], (E2 - mulst[l]) / tlst[l],
                                     dlst[l, 0] / tlst[l], dlst[l, 1] / tlst[l], sign=-1)
                    temp1 = PI * t2.real * f - t2.imag * phi0
                    temp2 = t2.real * phi0 + PI * t2.imag * f
                    kh.add_matrix_element(temp1, l, a2, a2, charge - 1, a1, b1, charge, 3)
                    kh.add_matrix_element(temp2, l, a2, a2, charge - 1, a1, b1, charge, 4)
        # Loop over final state, conserving charge
        for a2 in statesdm[charge]:
            if charge != ncharge - 1:
                # Intermediate state in higher charge state
                for c in statesdm[charge + 1]:
                    if b1 == a2:  # vertices on upper prop -> state on lower prop cannot change
                        E1 = E[c] - E[a2]
                        for l in range(nleads):
                            t2 = Tba[l, c, a1] * Tba[l, c, a2].conj()
                            f = fermi_func((E1 - mulst[l]) / tlst[l])
                            phi0 = phi((E1 - mulst[l]) / tlst[l], dlst[l, 0] / tlst[l],
                                       dlst[l, 1] / tlst[l], sign=1)
                            temp1 = -PI * t2.real * f - t2.imag * phi0
                            temp2 = -PI * t2.imag * f + t2.real * phi0
                            kh.add_matrix_element(temp1, l, a2, a2, charge, a1, b1, charge, 3)
                            kh.add_matrix_element(temp2, l, a2, a2, charge, a1, b1, charge, 4)
                    if a1 == a2:  # vertices on lower prop -> state on upper prop cannot change
                        E1 = E[c] - E[a2]
                        for l in range(si.nleads):
                            t2 = Tba[l, c, a2] * Tba[l, c, b1].conj()
                            f = fermi_func((E1 - mulst[l]) / tlst[l])
                            phi0 = phi((E1 - mulst[l]) / tlst[l], dlst[l, 0] / tlst[l],
                                       dlst[l, 1] / tlst[l], sign=1)
                            temp1 = - PI * t2.real * f + t2.imag * phi0
                            temp2 = - PI * t2.imag * f - t2.real * phi0
                            kh.add_matrix_element(temp1, l, a2, a2, charge, a1, b1, charge, 3)
                            kh.add_matrix_element(temp2, l, a2, a2, charge, a1, b1, charge, 4)
            if charge != 0:
                # Intermediate state in lower charge state
                for c in statesdm[charge - 1]:
                    if b1 == a2:  # vertices on upper prop -> state on lower prop cannot change
                        E1 = E[a2] - E[c]
                        for l in range(nleads):
                            t2 = Tba[l, a2, c] * Tba[l, a1, c].conj()
                            f = fermi_func(-(E1 - mulst[l]) / tlst[l])
                            phi0 = phi((E1 - mulst[l]) / tlst[l], dlst[l, 0] / tlst[l],
                                       dlst[l, 1] / tlst[l], sign=-1)
                            temp1 = - PI * t2.real * f + t2.imag * phi0
                            temp2 = - PI * t2.imag * f - t2.real * phi0
                            kh.add_matrix_element(temp1, l, a2, a2, charge, a1, b1, charge, 3)
                            kh.add_matrix_element(temp2, l, a2, a2, charge, a1, b1, charge, 4)
                    if a1 == a2:  # vertices on lower prop -> state on upper prop cannot change
                        E1 = E[a2] - E[c]
                        for l in range(nleads):
                            t2 = Tba[l, b1, c] * Tba[l, a2, c].conj()
                            f = fermi_func(-(E1 - mulst[l]) / tlst[l])
                            phi0 = phi((E1 - mulst[l]) / tlst[l], dlst[l, 0] / tlst[l],
                                       dlst[l, 1] / tlst[l], sign=-1)
                            temp1 = - PI * t2.real * f - t2.imag * phi0
                            temp2 = - PI * t2.imag * f + t2.real * phi0
                            kh.add_matrix_element(temp1, l, a2, a2, charge, a1, b1, charge, 3)
                            kh.add_matrix_element(temp2, l, a2, a2, charge, a1, b1, charge, 4)

    def generate_col_nondiag_kern_1st_order_nd(self, a1, charge):
        r""" Calculates a column in :math:`W_{nd}^{(1)}`, the part of the full first order off-diagonal kernel
        connecting non-diagonal :math:`|a2><b2|` and diagonal :math:`|a1><a1|` entries of the density matrix.
        Separates the real and imaginary parts into two matrices. Assumes that the wide band limit is valid.

        Parameters
        ----------
        a1 : int
            state indexing the column

        charge: int
            charge of state a1

        self.ImWdn : ndarray
            (Modifies) :math:`\Im(W_{nd}^{(1)})`, the imaginary part of the first order kernel that connects
            non-diagonal and diagonal elements of the density matrix.
        self.ReWdn : ndarray
            (Modifies) :math:`\Re(W_{nd}^{(1)})`, the real part of the first order kernel that connects
            non-diagonal and diagonal elements of the density matrix.
        """
        # Variables naming follow Leijnse et al.PRB 78, 235424(2008).

        (E, si, Tba, mulst, tlst, dlst) = (self.qd.Ea, self.si, self.leads.Tba,
                                           self.leads.mulst, self.leads.tlst, self.leads.dlst)

        si, kh = self.si, self.kernel_handler
        nleads, statesdm, ncharge = si.nleads, si.statesdm, si.ncharge
        PI = np.pi

        if charge != ncharge - 1:
            # Loop over final state, adding electron to the QD
            for a2 in statesdm[charge + 1]:
                E2 = E[a2] - E[a1]
                for b2 in statesdm[charge + 1]:
                    if a2 == b2:  # Final state must be off-diagonal
                        continue
                    E1 = E[b2] - E[a1]
                    for l in range(nleads):
                        t2 = Tba[l, a2, a1] * Tba[l, b2, a1].conj()
                        f = fermi_func((E1 - mulst[l]) / tlst[l]) + fermi_func((E2 - mulst[l]) / tlst[l])
                        phi0 = delta_phi((E1 - mulst[l]) / tlst[l], (E2 - mulst[l]) / tlst[l],
                                         dlst[l, 0] / tlst[l], dlst[l, 1] / tlst[l])
                        temp1 = PI * t2.real * f - t2.imag * phi0
                        temp2 = t2.real * phi0 + PI * t2.imag * f
                        kh.add_matrix_element(temp1, l, a2, b2, charge + 1, a1, a1, charge, 5)
                        kh.add_matrix_element(temp2, l, a2, b2, charge + 1, a1, a1, charge, 6)
        if charge != 0:
            # Loop over final states, removing electron from the QD
            for a2 in statesdm[charge - 1]:
                E1 = E[a1] - E[a2]
                for b2 in statesdm[charge - 1]:
                    if b2 == a2:  # Final state must be off-diagonal
                        continue
                    E2 = E[a1] - E[b2]
                    for l in range(si.nleads):
                        t2 = Tba[l, a1, b2] * Tba[l, a1, a2].conj()
                        f = fermi_func(-(E1 - mulst[l]) / tlst[l]) + fermi_func(-(E2 - mulst[l]) / tlst[l])
                        phi0 = delta_phi((E1 - mulst[l]) / tlst[l], (E2 - mulst[l]) / tlst[l],
                                         dlst[l, 0] / tlst[l], dlst[l, 1] / tlst[l], sign=-1)
                        temp1 = PI * t2.real * f - t2.imag * phi0
                        temp2 = t2.real * phi0 + PI * t2.imag * f
                        kh.add_matrix_element(temp1, l, a2, b2, charge - 1, a1, a1, charge, 5)
                        kh.add_matrix_element(temp2, l, a2, b2, charge - 1, a1, a1, charge, 6)
        # Loop over final state conserving charge
        for a2 in statesdm[charge]:
            for b2 in statesdm[charge]:
                if a2 == b2:  # Final state must be off-diagonal
                    continue
                if charge != ncharge - 1:
                    # Intermediate state in higher charge state
                    for c in statesdm[charge + 1]:
                        if a1 == b2:
                            E1 = E[c] - E[b2]
                            for l in range(nleads):
                                t2 = Tba[l, c, a1] * Tba[l, c, a2].conj()
                                f = fermi_func((E1 - mulst[l]) / tlst[l])
                                phi0 = phi((E1 - mulst[l]) / tlst[l], dlst[l, 0] / tlst[l], dlst[l, 1] / tlst[l])
                                temp1 = - PI * t2.real * f - t2.imag * phi0
                                temp2 = - PI * t2.imag * f + t2.real * phi0
                                kh.add_matrix_element(temp1, l, a2, b2, charge, a1, a1, charge, 5)
                                kh.add_matrix_element(temp2, l, a2, b2, charge, a1, a1, charge, 6)
                        if a1 == a2:
                            E1 = E[c] - E[a2]
                            for l in range(nleads):
                                t2 = Tba[l, c, b2] * Tba[l, c, a1].conj()
                                f = fermi_func((E1 - mulst[l]) / tlst[l])
                                phi0 = phi((E1 - mulst[l]) / tlst[l], dlst[l, 0] / tlst[l], dlst[l, 1] / tlst[l])
                                temp1 = - PI * t2.real * f + t2.imag * phi0
                                temp2 = - PI * t2.imag * f - t2.real * phi0
                                kh.add_matrix_element(temp1, l, a2, b2, charge, a1, a1, charge, 5)
                                kh.add_matrix_element(temp2, l, a2, b2, charge, a1, a1, charge, 6)
                if charge != 0:
                    # Intermediate state in lower charge state
                    for c in statesdm[charge - 1]:
                        if a1 == b2:
                            E1 = E[b2] - E[c]
                            for l in range(nleads):
                                t2 = Tba[l, a2, c] * Tba[l, a1, c].conj()
                                f = fermi_func(-(E1 - mulst[l]) / tlst[l])
                                phi0 = phi((E1 - mulst[l])/tlst[l], dlst[l, 0] / tlst[l], dlst[l, 1] / tlst[l], sign=-1)
                                temp1 = - PI * t2.real * f + t2.imag * phi0
                                temp2 = - PI * t2.imag * f - t2.real * phi0
                                kh.add_matrix_element(temp1, l, a2, b2, charge, a1, a1, charge, 5)
                                kh.add_matrix_element(temp2, l, a2, b2, charge, a1, a1, charge, 6)
                        if a1 == a2:
                            E1 = E[a2] - E[c]
                            for l in range(nleads):
                                t2 = Tba[l, a1, c] * Tba[l, b2, c].conj()
                                f = fermi_func(-(E1 - mulst[l]) / tlst[l])
                                phi0 = phi((E1 - mulst[l])/tlst[l], dlst[l, 0] / tlst[l], dlst[l, 1] / tlst[l], sign=-1)
                                temp1 = - PI * t2.real * f - t2.imag * phi0
                                temp2 = - PI * t2.imag * f + t2.real * phi0
                                kh.add_matrix_element(temp1, l, a2, b2, charge, a1, a1, charge, 5)
                                kh.add_matrix_element(temp2, l, a2, b2, charge, a1, a1, charge, 6)

    def generate_row_inverse_Liouvillian(self, a1, b1, charge):
        """ Calculates a row of :math:`1/(L_{nn})` (in practice only the diagonal is needed) where :math:`L_{nn}` is the part
        of the Liouvillian connecting non-diagonal states :math:`|i><j|` with non-diagonal states :math:`|k><l|`.

        Parameters
        ----------
        a1 : int
            first state indexing the column

        b1 : int
            second state indexing the column

        charge: int
            charge of state a1 and b1

        self.Wnn : ndarray
            (Modifies) an ndarray representing :math:`1/L_{nn}`. This kernel has noffdiag * noffdiag entries.

        """
        # Variables names follow Leijnse et al PRB 78, 235424 (2008).
        E = self.qd.Ea

        minE = 1e-10
        E1 = E[a1] - E[b1]
        if minE > E1 >= 0:
            E1 = minE
        elif -minE < E1 <= 0:
            E1 = -minE

        self.kernel_handler.add_element_Lnn(a1, b1, charge, 1.0/E1)

    def set_Ozaki_params(self):
        """
        Generates and stores the residues and poles for the Ozaki expansion
        of tanh(z).

        Parameters
        ----------
        self.BW_Ozaki_expansion : double
            (Modifies) The band width (over temperature) for which the Ozaki poles and residues are generated.

        self.Ozaki_poles_and_residues : ndarray
            a ndarray where the first column contains the reciprocal of the poles
            and the second columns contains the residues.
        """
        BW_T = (abs(self.leads.dlst[0][0]) + abs(self.leads.dlst[0][1])) / 2.0 / min(self.leads.tlst)
        if self.BW_Ozaki_expansion < BW_T:
            self.Ozaki_poles_and_residues = BW_Ozaki(BW_T)
            self.BW_Ozaki_expansion = BW_T

    def generate_vec(self, phi0):
        """"""
        raise NotImplementedError('Matrix free methods are not supported by the RTD approach.')
import copy
import numpy as np
from ..wrappers.mytypes import longnp
from ..approach.aprclass import Approach
from ..indexing import StateIndexingDM
from ..indexing import StateIndexingDMc
from ..qdot import QuantumDot
from ..leadstun import LeadsTunneling
from .funcprop import FunctionProperties
from .various import get_phi0
from .various import get_phi1
from .various import print_state
from .various import print_all_states
from .various import sort_eigenstates
from .various import remove_coherences
from .various import remove_states
from .various import use_all_states
from .validation import validate_kerntype
from .validation import validate_itype
from .validation import validate_indexing
from ..approach.base.pauli import ApproachPauli as ApproachPyPauli
from ..approach.base.lindblad import ApproachLindblad as ApproachPyLindblad
from ..approach.base.redfield import ApproachRedfield as ApproachPyRedfield
from ..approach.base.neumann1 import Approach1vN as ApproachPy1vN
from ..approach.base.neumann2 import Approach2vN as ApproachPy2vN
from ..approach.base.RTD import ApproachPyRTD as ApproachPyRTD
    from ..approach.base.c_pauli import ApproachPauli
    from ..approach.base.c_lindblad import ApproachLindblad
    from ..approach.base.c_redfield import ApproachRedfield
    from ..approach.base.c_neumann1 import Approach1vN
    from ..approach.base.c_neumann2 import Approach2vN
    from ..approach.base.c_RTD import ApproachRTD
            sub_class = getattr(self, attribute_map[item])
            setattr(sub_class, item, value)

    def change_si(self):
        si = self.si
        icn = self.Approach.indexing_class_name
        # noinspection PyTypeHints
        if not isinstance(si, self.globals[icn]):
            self.si = self.globals[icn](si.nsingle, si.indexing, si.symmetry, si.nleads)
            self.qd.si = self.si
            self.leads.si = self.si
        self.si.states_changed = True

    # kerntype
    def get_kerntype(self):
        return self.appr.kerntype

    def set_kerntype(self, value):
        if isinstance(value, str):
            if self.appr.kerntype != value:
                approach_string = value[0].capitalize() + value[1:]
                self.Approach = self.globals['Approach'+approach_string]
                self.change_si()
                self.appr = self.Approach(self)
        else:
            if issubclass(value, Approach):
                self.Approach = value
                self.change_si()
                self.appr = self.Approach(self)
    kerntype = property(get_kerntype, set_kerntype)

    # indexing
    def get_indexing(self):
        return self.si.indexing

    def set_indexing(self, value):
        if self.si.indexing != value:
            print('WARNING: Cannot change indexing from \''+self.si.indexing+'\' to \''+value
                  + '\' Consider constructing a new system using \'indexing='+value+'\'')
    indexing = property(get_indexing, set_indexing)

    # dband
    def get_dband(self):
        return self.funcp.dband

    def set_dband(self, value):
        self.funcp.dband = value
        self.leads.change(dlst=value)
    dband = property(get_dband, set_dband)

    def add(self, hsingle=None, coulomb=None, tleads=None, mulst=None, tlst=None, dlst=None):
        """
        Adds the values to the specified dictionaries and correspondingly redefines
        relevant many-body properties of the system.

        Parameters
        ----------
        hsingle,coulomb,tleads,mulst,tlst,dlst : dict
            Dictionaries describing what values to add.
            For example, tleads[(lead, state)] = value to add.
        """
        if not (hsingle is None and coulomb is None):
            self.qd.add(hsingle, coulomb)
        if not (tleads is None and mulst is None and tlst is None and dlst is None):
            self.leads.add(tleads, mulst, tlst, dlst)

    def change(self, hsingle=None, coulomb=None, tleads=None, mulst=None, tlst=None, dlst=None):
        """
        Changes the values of the specified dictionaries and correspondingly redefines
        relevant many-body properties of the system.

        Parameters
        ----------
        hsingle,coulomb,tleads,mulst,tlst,dlst : dict
            Dictionaries describing what values to change.
            For example, tleads[(lead, state)] = value to change.
        """
        if not (hsingle is None and coulomb is None):
            self.qd.change(hsingle, coulomb)
        if not (tleads is None and mulst is None and tlst is None and dlst is None):
            self.leads.change(tleads, mulst, tlst, dlst)

    def get_phi0(self, b, bp):
        """
        Get the reduced density matrix element corresponding to
        many-body states b and bp.
        """
        return get_phi0(self, b, bp)

    def get_phi1(self, l, c, b):
        """
        Get the energy integrated current amplitudes corresponding to
        lead l and many-body states c and b.
        """
        return get_phi1(self, l, c, b)

    def print_state(self, b, eps=0.0, prntq=True, filename=None, separator=''):
        """
        Prints properties of given many-body eigenstate of the quantum dot Hamiltonian
        """
        print_state(self, b, eps, prntq, filename, separator)

    def print_all_states(self, filename, eps=0.0, separator='', mode='w'):
        """
        Prints properties of all many-body eigenstates to a file.
        """
        print_all_states(self, filename, eps, separator, mode)

    def sort_eigenstates(self, srt=None):
        """
        Sort many-body states of the system by given order of properties.
        """
        sort_eigenstates(self, srt=srt)

    def remove_coherences(self, dE):
        """
        Remove the coherences with energy difference larger than dE.
        """
        remove_coherences(self, dE)


    def remove_states(self, dE):