Compact algebraic equations defining the equation-of-motion coupled-cluster (EOM-CC) methods for ionization potentials (IP-EOM-CC) have been derived and computer implemented by virtue of a symbolic algebra system largely automating these processes. Models with connected cluster excitation operators truncated after double, triple, or quadruple level and with linear ionization operators truncated after two-hole-one-particle (2h1p), three-hole-two-particle (3h2p), or four-hole-three-particle (4h3p) level (abbreviated as IP-EOM-CCSD, CCSDT, and CCSDTQ, respectively) have been realized into parallel algorithms taking advantage of spin, spatial, and permutation symmetries with optimal size dependence of the computational costs. They are based on spin-orbital formalisms and can describe both alpha and beta ionizations from open-shell (doublet, triplet, etc.) reference states into ionized states with various spin magnetic quantum numbers. The application of these methods to Koopmans and satellite ionizations of N-2 and CO (with the ambiguity due to finite basis sets eliminated by extrapolation) has shown that IP-EOM-CCSD frequently accounts for orbital relaxation inadequately and displays errors exceeding a couple of eV. However, these errors can be systematically reduced to tenths or even hundredths of an eV by IP-EOM-CCSDT or CCSDTQ. Comparison of spectroscopic parameters of the FH+ and NH+ radicals between IP-EOM-CC and experiments has also underscored the importance of higher-order IP-EOM-CC treatments. For instance, the harmonic frequencies of the A (2)Sigma(-) state of NH+ are predicted to be 1285, 1723, and 1705 cm(-1) by IP-EOM-CCSD, CCSDT, and CCSDTQ, respectively, as compared to the observed value of 1707 cm(-1). The small adiabatic energy separation (observed 0.04 eV) between the X (2)Pi and a (4)Sigma(-) states of NH+ also requires IP-EOM-CCSDTQ for a quantitative prediction (0.06 eV) when the a (4)Sigma(-) state has the low-spin magnetic quantum number (s(z)=1/2). When the state with s(z)=3/2 is sought, the energy separations converge much more rapidly with the IP-EOM-CCSD value (0.03 eV) already being close to the observed (0.04 eV). (c) 2006 American Institute of Physics.