Below I explain how to enumerate all decks without bad pairs explicitly. The exact number turns out to be $3668033946384704437729512814619767610579526911188666362431432294400.$ Dividing by $52!$, I get roughly $0.045476282331.$
We are going to write a recurrence relation for $f(H,k)$, which is the number of partial decks with $t$ cards appearing $H_t$ times, the last card making its $k$th appearance. The total number of decks is then $f((0,0,0,0,13),4)$.
The base case is $f((12,1,0,0,0),1) = 52$. Now suppose we're given some $H$ and $k$ such that $H_k \geq 1$. First form H' by moving one "rank" from $k$ to $k-1$. Now for each $t \neq k - 1$, we can extend each of the $f(H_k,t)$ partial decks to a new deck in $H_k(k-1) \cdot (4-(k-1))$ ways. For $t = k-1$, we can extend the $f(H_k,k-1)$ partial decks in only $(H_k(k-1)-1) \cdot (4-(k-1))$ ways (since one rank is taken).
There are only 9481 pairs $H,k$ to consider, so even my lousy python implementation which uses a hash table to keep track of all those pairs is fast enough.
Michael Lugo's method gives $(1-3/51)^{51} \approx 0.0454176$, which is quite close to the right answer. The correct probability is slightly larger since if we imagine a process in which cards are revealed one by one, the probability that the next card is bad starts at $3/51$ but eventually decreases (roughly).
His approximation $e^{-3} \approx 0.049787$ is much too big, and that's because $(1-3/n)^n$ approaches $e^{-3}$ from below.
leonbloy's approximation is also too small, for similar reasons. We can imagine a process in which ranks are revealed one by one. After the first rank is revealed, a smaller number of positions are adjacent, and so the probability decreases.