The 6-sided dice case can be encoded as a Markov chain on the state space $\{0,1,2,3,\partial\}$ where each state $0\leqslant i\leqslant3$ means that $i$ pairs from the set $\{\{1,6\},\{2,5\},\{3,4\}\}$ have already been visited but none is completed yet, and state $\partial$ means that at least one of these pairs has been completed. For example, after the results 26622 the state of the chain is $2$, after the results 266224246 the state of the chain is $3$, and after the results 2662242465 the state of the chain is $\partial$.
The probability that any two dice add to 7 on a throw of $n$ is $P_0[T\leqslant n]$ where $T$ is the hitting time of $\partial$ and the subscript $0$ refers to the starting point of the chain.
To characterize $P_0[T\leqslant n]$ for every $n$, one usually computes the generating functions $u_i=E_i[s^T]$ for $0\leqslant i\leqslant3$, where $|s|\leqslant1$. Writing down carefully the transition probabilities of the chain, one sees that the Markov property after one step yields the relations
$$
u_0=su_1,\quad u_1=s(\tfrac16u_1+\tfrac46u_2+\tfrac16),\quad u_2=s(\tfrac26u_2+\tfrac26u_3+\tfrac26),\quad u_3=s(\tfrac36u_3+\tfrac36).
$$
Solving this Cramér system, one gets (something similar to)
$$
u_0=\frac{s^2(6+3s+s^2)}{(2-s)(3-s)(6-s)}.
$$
Thus, $u_0$ is a rational fraction with respect to $s$, which can be decomposed as
$$
u_0=\frac{c_2}{2-s}+\frac{c_3}{3-s}+\frac{c_6}{6-s}+as+b,
$$
for some suitable constants $a$, $b$, $c_2$, $c_3$ and $c_6$. Expanding each fraction as a power series in $s$ and collecting all the terms of degree at least $n$ yields, for every $n\geqslant2$,
$$
P[T\geqslant n]=\sum_k\frac{c_k}{k-1}\frac1{k^n},
$$
where the sum runs over $k$ in $\{2,3,6\}$.
If one uses some $2z$-sided dice instead, one finds similarly
$$
P[T\geqslant n]=\sum_k\frac{c_k}{k-1}\frac1{k^n},
$$
for some suitable constants $c_k$, where the sum runs over $k=2z/\ell$ with $1\leqslant\ell\leqslant z$. In particular, when $n\to\infty$,
$$
P[T\geqslant n]\sim\frac{c_2}{2^n}.
$$
Going back to the 6-sided case, $c_2=\lim\limits_{s\to2}(2-s)u_0=16$ hence
$$
P[T\geqslant n]\sim\frac{16}{2^{n}}.
$$
Best Answer
Concerning a "closed" ( = finite summation) formula, start from $$ \eqalign{ & N_b (s - n,m - 1,n) = {\rm No}{\rm .}\,{\rm of}\,{\rm solutions}\,{\rm to}\;\left\{ \matrix{ {\rm 1} \le {\rm integer}\;x_{\,j} \le m \hfill \cr x_{\,1} + x_{\,2} + \; \cdots \; + x_{\,n} = s \hfill \cr} \right. = \cr & = {\rm No}{\rm .}\,{\rm of}\,{\rm solutions}\,{\rm to}\;\left\{ \matrix{ 0 \le {\rm integer}\;y_{\,j} \le m - 1 \hfill \cr y_{\,1} + y_{\,2} + \; \cdots \; + y_{\,n} = s - n \hfill \cr} \right. \cr} $$ where $N_b$ is given by $$ N_b (s,r,m)\quad \left| {\;0 \leqslant \text{integers }s,m,r} \right.\quad = \sum\limits_{\left( {0\, \leqslant } \right)\,\,k\,\,\left( { \leqslant \,\frac{s}{r+1}\, \leqslant \,m} \right)} {\left( { - 1} \right)^k \binom{m}{k} \binom { s + m - 1 - k\left( {r + 1} \right) } { s - k\left( {r + 1} \right)}\ } $$ as explained in this related post.
Then the number of ways to obtain a sum $x \le s$ is given by $$ \eqalign{ & M(x,m,n) = \sum\limits_{x\, \le \,s\,\left( { \le \,m\,n} \right)\,} {N_b (s - n,m - 1,n)} = \cr & = m^{\,n} - \sum\limits_{0\, \le \,s\, \le \,x - 1\,} {N_b (s - n,m - 1,n)} = \cr & = m^{\,n} - \sum\limits_{0\, \le \,s\, \le \,x - n - 1\,} {N_b (s,m - 1,n)} = \cr & = m^{\,n} - \sum\limits_{0\, \le \,s\, \le \,x - n - 1\,} {\sum\limits_{\left( {0\, \le } \right)\,\,k\,\,\left( { \le \,{s \over m}\, \le \,n} \right)} {\left( { - 1} \right)^k \left( \matrix{ n \hfill \cr k \hfill \cr} \right)\left( \matrix{ s + n - 1 - k\,m \cr s - k\,m \cr} \right)} } = \cr & = m^{\,n} - \sum\limits_{\left( {0\, \le } \right)\,\,k\,\,\left( { \le \,{s \over m}\, \le \,n} \right)} {\left( { - 1} \right)^k \left( \matrix{ n \hfill \cr k \hfill \cr} \right)\left( \matrix{ x - 1 - k\,m \cr x - n - 1 - k\,m \cr} \right)} \cr} $$
and in fact $M(90,20,5)=3003$.
Note that, as explained in the cited related post, the problem has the geometric equivalent of finding:
the number of integral points on the diagonal plane $y_1, \cdots, y_n=s-n$, intercepted by a $n$-D cube with side $[0,m-1]$.
The formula for $N_b$ corresponds to calculating the points on the whole plane ($k=0$) and subtracting those pertaining to the surrounding cubes.
The geometric analogy clearly shows that $N_b(nm-s,m,n)=N_b(s,m,n)$.