A method is developed to compute the nonadiabatic coupling vectors (NACVs) between electronic ground and excited states as well as between any possible pair of excited states within the framework of linear response time-dependent density functional theory (TDDFT) in the adiabatic approximation. The development is an extension to our previous work on surface hopping dynamics [E. Tapavicza , Phys. Rev. Lett. 98, 023001 (2007)] for which we improve the description of the TDDFT approximation of the excited state wavefunctions by means of linear response orbitals. The method is first validated on the H+H-2 system that has a region of strong coupling near the conical intersection at the equilateral geometry. These results confirm the quality and the numerical efficiency of the approach, which has an accuracy comparable to the one achieved with wavefunction-based methods. Finally, we apply the method to the calculation of the NACVs of protonated formaldimine (NH2CH2+) along a surface hopping trajectory initiated in the second excited state.