  We report a calculation of the perturbative matching coefficients for the
transverse-momentum-dependent parton distribution functions for quark at the
next-to-next-to-next-to-leading order in QCD, which involves calculation of
non-standard Feynman integrals with rapidity divergence. We introduce a set of
generalized Integration-By-Parts equations, which allows an algorithmic
evaluation of such integrals using the machinery of modern Feynman integral
calculation.
