Probability for hand with no sets (3+ of a kind / straight flush of 3+ cards) or potential sets (no additional card can make a set) in a card game? Problem
In a gin variant, you win immediately if you are dealt a hand with no sets (3+ of a kind or suited run of 3+ cards) and no "potential sets", i.e. no additional card added to the hand could make a set. What's the probability of this happening for a hand of $n$ cards? I actually play this game ($n=9$) and my friend and I were curious about this probability, which is the motivation for this question. It's easy to compute numerically for any $n$ and some special $n$ values are also easy (e.g. $n$=1, $n>13$) but an analytic solution for general $n$ seems deceptively difficult.
(Bad) Approach Idea
I cannot seem to find a nice way to do it analytically at all. Obviously we want the total number of ways to make this special hand, $N_h$, divided by $N={52 \choose n}$ possible hands. 
The no 3+-of-a-kind/potential 3+-of-a-kind requirement is simple enough because a potential 3-of-a-kind is just a pair. This means we must select $n$ distinct values from the $13$ values in the deck $\to {13 \choose n}$.
However the suited runs are much more difficult. Depending on the distribution of the $n$ values across the $4$ suits, the number of ways of making runs/potential runs vary. The only thing I can think of is to enumerate the probabilities of runs/potential runs for each number of cards in a suit. E.g. $1$ card in a suit never makes a run/potential run. $2$ cards has $n-m+1$, $m$ card runs and $n-2$ "gapped" potential runs. $n=3$ then is more complicated with the placement of the $3$rd relative to the other $2$. There are also complications with double counting since these things are not mutually exclusive. I have some vague notions that these enumerations have some relationship to the polyhedral numbers but haven't pursued that in depth because I'm convinced there's a better approach and I'm hoping someone here can help me find it!
Numerics
It's easy to write a program to do this numerically for some $n$. Just check that there are no pairs and that within a suit all cards are at least $2$ away from every other card. Here's a quick Mathematica code that does it
deal[n_] := QuotientRemainder[RandomSample[Range[0, 51], n], 13]
valid[hand_] := DuplicateFreeQ[Transpose[hand][[-1]]] &&
    Min[GroupBy[hand,First,Min[Abs[Differences/@Subsets[Last/@#,{2}]]]&]]>2
prob[n_, num_] := N@Count[Table[valid[deal[n]], {i, num}], True]/num

For the game I play $n=9$ and dealing half a billion hands gives $P\approx 0.00185276$
 A: Consider the sequence of $n$ distinct ranks in a hand without potential suited runs. Let $j$ be the number of groups of ranks separated by gaps of length at least $2$. Let $k$ be the number of contiguous runs of ranks. And let $m$ be the number of contiguous runs of ranks of length at least $2$.
Consider the number of suit options for each rank, starting from the highest. The first rank can have all $4$ suits, and any rank following a gap of length at least $2$ can also have all $4$ suits, so there are $j$ ranks that can have all $4$ suits. Other ranks can usually have one of $3$ suits, unless they’re in the third or later position of a run of contiguous ranks, in which case they can have only one of $2$ suits. So there are $j$ ranks with $4$ options, $k+m-j$ ranks with $3$ options and $n-(k+m)$ ranks with $2$ options. The positions of the $m$ runs of length at least $2$ can be chosen among the $k$ runs in $\binom km$ ways. The $n-(k+m)$ extra ranks (the ones with $2$ options) can be distributed over the $m$ runs of length at least $2$ in $\binom{n-(k+m)+m-1}{m-1}=\binom{n-k-1}{m-1}$ ways. There are $\binom{k-1}{j-1}$ ways to choose the $j-1$ long gaps among the $k-1$ gaps. And the $13-n-(k-1)-(j-1)$ ranks that are not selected and not yet accounted for in the gaps can be distributed over $j+1$ slots in $\binom{13-n-(k-1)-(j-1)+j}j=\binom{15-n-k}j$ ways. Thus in total there are
$$
\sum_{k=1}^{\min(n,15-n)}\sum_{m=0}^{n-k}\sum_{j=1}^k\binom km\binom{n-k-1}{m-1}\binom{k-1}{j-1}\binom{15-n-k}j4^j3^{k+m-j}2^{n-(k+m)}
$$
hands without potential suited runs. Here’s code to evaluate this in Sage:
j,k,m,n=var('j,k,m,n')
def binom(n,k):
    return 1 if n == -1 and k == -1 else binomial(n,k);
sum(binomial (k,m) * binom (n-k-1,m-1) * binomial (k-1,j-1) * binomial (15-n-k,j) * 4**j * 3**(k+m-j) * 2**(n-(k+m)) for k in range(1,min(n+1,16-n)) for j in range (1,n+1) for m in range(0,n-k+1))

(Note that the case $m=0$ and thus $k=n$ is treated correctly, since for this case $\binom{-1}{-1}=1$, but this is not correctly implemented in Sage; hence the custom function binom.)
For $n=9$ the result is $6820068$, so the probability to draw such a hand is
$$
\frac{6820068}{\binom{52}9}=\frac{1705017}{919768850}\approx0.0018537\;,
$$
in good agreement with your numerical result (which you stated with excess precision, as its standard deviation is $\sqrt{\frac{p(1-p)}n}\approx2\cdot10^{-6}$).
