Revised
Interesting question.
Here's a thought:
You can think of a ring, such as $\mathbb Z$, in terms of its monoid of affine endomorphisms
$x \rightarrow a x + b$. The action of this monoid, together with a choice for 0 and 1, give the structure of the ring. However, the monoid is not finitely generated, since the
multiplicative monoid of $\mathbb Z$ is the free abelian monoid on the the primes, times
the order 2 group generated by $-1$.
If you take a submonoid that uses only
one prime, it is quasi-isometric to a quotient of the hyperbolic plane by an action of $\mathbb Z$, which is multiplication by $p$ in the upper half-space model. To see this, place a dots labeled by integer $n$ at position $(n*p^k, p^k)$ in the upper half plane, for every pair of integers $(n,k)$, and connect them by horizontal line segments and by vertical line segments whenever points are in vertical alignment. The quotient of upper half plane by the hyperbolic isometry $(x,y) \rightarrow (p*x, p*y)$ has a copy of the Cayley graph for this monoid. This is also quasi-isometric to the 1-point union of two copies of the hyperbolic plane, one for negative integers, one for positive integes. It's a fun exercise,
using say $p = 2$. Start from 0, and recursively build the graph by connecting $n$ to $n+1$ by one color arrow, and $n$ to $2*n$ by another color arrow. If you arrange positive integers in a spiral, you can make a neat drawing of this graph (or the corresponding graph for a different prime.) The negative integers look just the same, but with the successor arrow reversed.
If you use several primes, the picture gets more complicated. In any case, one can take rescaled limits of these graphs, based at sequence of points, and get asymptotic cones for the monoid. The graph is not homogeneous, so there is not just one limit.
Another point of view is to take limits of $\mathbb Z$ without rescaling, but
with a $k$-tuple of constants $(n_1, \dots , n_k)$. The set of possible identities among polynomials in $k$ variables is compact, so there is a compact space of limit rings for $\mathbb Z$ with $k$ constants. Perhaps this is begging the question: the identitites that define the limits correspond to diophantine equations that have infinitely many solutions.
Rescaling may eliminate some of this complexity.
A homomorphism $\mathbb Z[x,y,\dots,z]/P$ to
$\mathbb Z$ gives a homomomorphism of the corresponding monoids, so an infinite sequence of these gives an action on some asymptotic cone for the affine monoid for $\mathbb Z$.
With the infinite set of primes, there are other plausible choices for how to define length; what's the best choice depends on whether and how one can prove anything of interest.
The number of rational solutions to your equation is finite. In short: your equation defines a genus $3$ curve, as follows from a straightforward computation and an application of Riemann--Hurwitz; finally, by Faltings' theorem, the number of rational points on a curve of genus $>1$ is finite.
One shows this as follows. Your equation defines a smooth projective curve $C$ whose function field $K$ is generated by $p$ and $q$. Take the second equation
$$
\frac{1}{\sqrt{p-1}} - \frac{1}{\sqrt{p+1}} = \sqrt{\frac{q}{q^2-1}}.
$$
Squaring it, we get
$$
\frac{q}{q^2-1} = \frac{2p}{p^2-1}-\frac{2}{\sqrt{p^2-1}}
$$
showing that $r:=\sqrt{p^2-1}$ is in $K$. Put differently, $C$ maps to the rational curve $C_0$ given by
$$
p^2-r^2=1.
$$
In fact, since as we've seen $q/(q^2-1) = f$ with $f=2p/r^2-2/r$, the curve $C$ is a double cover of $C_0$ (and therefore hyperelliptic). We determine the number of branch points.
Rewriting the last equation, we find
$$
( qf ) ^2- (qf) - f^2 = 0,
$$
which is a quadratic (in variable $qf$) with discriminant $1+4f^2$, so $C \to C_0$ is ramified over the points on $C_0$ where $f = \pm i/2$. This last equality expands to $4(p-r)=\pm ir^2$, which gives $p+r=1/(p-r)=\mp 4i/r^2$, so $r=[(p+r)-(p-r)]/2=\mp (2i/r^2 + ir^2/8),$ giving four solutions for $r$ for each choice of sign, each corresponding to a unique value of $p$, so eight branch points $(p,r)$ in total. (Indeed, the eight values of $r$ are given by the zeros of $r^8 - 64r^6 + 32r^4 + 256
$, which sage factors as $(r^4 - 8r^3 + 16)(r^4 + 8r^3 + 16)$, so we indeed have $8$ distinct solutions.) Therefore $C$ has genus $\left\lfloor(8-1)/2\right\rfloor=3$ by Riemann--Hurwitz, which proves that your equation has finitely many solutions by Faltings' theorem.
Best Answer
I echo Richard Stanley in his pessimism for there being a nice formula or method for giving you the solutions. On the other hand, most of the solutions will have a_i being 0 for most of the coefficients i. Here is another way to look at it, which might help you write a program for computer search for small n.
Your constraint summing to $n$ gives that $a_i + \epsilon_i = \frac{n}{i}$, where for most $i$, $0<\epsilon_i < 2$. Since you know most $a_i$ will be zero, $s$ will be not far from a sum of terms of the form $\frac{n(i-1)}{2}$, so many solutions will be close to partitioning an integer near $\frac{2s}{n}$ into distinct parts (so e.g. a partition of 20 into {1,4,4,5,6} is not allowed). This roughening of the problem can be easily programmed without expanding the search space too much.
If you still hope for more, you might try deriving some recurrences which might help in a partial characterization of solutions. It is my feeling that a quick program can be developed that will give you the solution set for any feasible s even for n as large as 50.
Gerhard "Email Me About System Design" Paseman, 2011.07.11