We present recent developments in numerical methods for the energy density functional calculations, suitable for massively parallel computing: The green's function with shifted Krylov subspace method and the Fermi operator expansion method. Some applications of those methods to inhomogeneous nuclear matter are presented. This may elucidate the transport properties of the superfluid neutrons in the inner crust of neutron stars, which are essential ingredients in pulsar glitch phenomena.