According to http://www-math.mit.edu/~rstan/transparencies/snf.pdf,
the Smith Normal form for a nonsingular matrix $M$ is
$$
\textrm{Diag}(e_1,\ldots, e_n),
$$
where $e_1\cdots e_i = \gcd \big( i\times i\textrm{ minors of }M \big)$.

There are infinitely many $a_1,\ldots, a_n$ such that $M$ is nonsingular, so we assume that $M$ is nonsingular for convenience. The singular case may be treated separately.

**Case 1 : $M$ is nonsingular**

Since the matrix has $-1$ in most of the entries, the $i\times i$ minors are relatively easy to find.

From this, the $1\times 1$ entry $e_1$ is easily seen to be $1$.

The $2\times 2$ entry $e_2$ can be obtained from
$$
e_2=e_1e_2=\gcd\big( \{a_i+1 | i=1,\ldots,n\}, \{a_ia_j-1 | 1\leq i<j\leq n\}\big).
$$
Since $a_ia_j-1=(a_i+1)(a_j+1)-(a_i+1)-(a_j+1)$, we have
$$
e_2=\gcd\big( \{a_i+1 | i=1,\ldots,n\}\big).
$$
The $3\times 3$ entry $e_3$ is obtained from
$$
\begin{align}
e_1e_2e_3&=\gcd\big( \{(a_i+1)(a_j+1) | 1\leq i<j\leq n\}, \{(a_i+1)(a_j+1)(a_k+1) -(a_i+1)(a_j+1) - (a_j+1)(a_i+1) - (a_k+1)(a_i+1)|1\leq i<j<k\leq n\}\big)\\
&=\gcd\big( \{(a_i+1)(a_j+1) | 1\leq i<j\leq n\} \big)
\end{align}
$$

In general, for $2\leq i\leq n-1$,
$$
e_1\cdots e_i = \gcd\big( \{\prod_{k=1}^{i-1} (a_{j_k}+1) | 1\leq j_1< \cdots <j_{i-1}\leq n\}\big).
$$
Then $e_n$ is obtained from
$$
e_n=\frac{\det M}{e_1\cdots e_{n-1}}.$$

**Case 2: $M$ is singular**

This case, we have
$$
1 = \sum_{i=1}^n \frac1{a_i+1}.
$$
Since $1\neq \sum_{i=1}^{n-1}\frac1{a_i+1}$, the rank of $M$ is $n-1$ and we can proceed the same method as in Case 1, after putting all zeros on the last column and the last row by row/column operations.