First, let us compute $f(m)$, the number of ways to deal two cards to each of $m$ people which are all pairs. We can do this using a computer as follows. If we further let $f(m,r)$ be the probability that all $m$ people have pairs when dealt from a deck of $4r$ cards, then
$$
f(m,r) = r\binom42\big(f(m-1,r-1)+(m-1)f(m-2,r-1)\big)
$$
The first summand accounts for ways where the first player gets a pair that no one else has, and the second accounts for ways where some other player also gets that pair. The above recursive equation allows you to compute $f(m)=f(m,13)$ quickly using dynamic programming.
Next, use the generalized principle of inclusion exclusion. Letting $E_i$ be the event that the $i^{th}$ player has a pair, then
\begin{align}
P(n,k)
&=\sum_{i=0}^n(-1)^{i-k}\binom{i}{k}\binom{n}{i}\mathbb P(E_1\cap E_2\cap \dots \cap E_i)\\
&=\sum_{i=0}^n(-1)^{i-k}\binom{i}{k}\binom{n}{i}\frac{f(i)}{\frac{52!}{2^i(52-2i)!}}
\end{align}
You can see this in action here, just click the run button at the top and enter your desired value of $n$.
Also, there is the following "closed form" for $P(n,k)$:
$$
\boxed{P(n,k) = \binom{n}{k}\sum_{i=k}^n(-1)^{i-k}\binom{n-k}{i-k}\frac{2^i(52-2i)!}{52!}\cdot i!\cdot [x^i](1+6x+3x^2)^{13}}
$$
Here, the notation $[x^i]f(x)$ means the coefficient of $x^i$ in the polynomial $f(x)$.