The complete linearized Fokker–Planck collision operator has been implemented in the drift-kinetic code NEO (Belli and Candy 2008 Plasma Phys. Control. Fusion50 095010) for the calculation of neoclassical transport coefficients and flows. A key aspect of this work is the development of a fast numerical algorithm for treatment of the field particle operator. This Eulerian algorithm can accurately treat the disparate velocity scales that arise in the case of multi-species plasmas. Specifically, a Legendre series expansion in ξ (the cosine of the pitch angle) is combined with a novel Laguerre spectral method in energy to ameliorate the rapid numerical precision loss that occurs for traditional Laguerre spectral methods. We demonstrate the superiority of this approach to alternative spectral and finite-element schemes. The physical accuracy and limitations of more commonly used model collision operators, such as the Connor and Hirshman–Sigmar operators, are studied, and the effects on neoclassical impurity poloidal flows and neoclassical transport for experimental parameters are explored.