The solutions of Lyapunov and Sylvester equations can be represented as a sum of Hermitian matrices or sub-Gramians corresponding either to individual eigenvalues of the system or their pairwise combinations. The sub-Gramians corresponding to a given part of the spectrum of dynamical system can be used to estimate the modal interaction in this part of the spectrum. In this paper we propose an algorithm for the fast calculation of sub-Gramians for large electric power systems. Our numerical tests demonstrate that the sub-Gramians corresponding to the selected part of the spectrum can be calculated much faster than all sub-Gramians of the system. Therefore, the proposed algorithm can serve as a basis for the fast calculation of modal interactions in large-scale power systems in real time.