It can be done for the metric 

$$d(a,b)^2 = 1 -  \frac{(a,b)}{\sqrt{ab}},$$

and other similar ones like $d(a,b)^2 = 1 -  \frac{(a,b)^2}{ab}$, with some twists in the construction.

Suppose we want to embed $1,2,..., n$ in $\mathbb{R}^n$. We will first embed these in $\mathbb{R}^m$, where $m = lcm(1,2,...,n)$.

For each natural number $k\in\{1,...,n\}$ map it to the vector $v_k \in \mathbb{R}^m$ whose $i$-th entry is equal to $\sqrt{k}$ if $i$ and is a multiple of $k$ and $0$ otherwise. Noticing that the vectors $v_a$ and $v_b$ are only both non-zero at the entries multiple of $[a,b] = lcm(a,b)$ we get:

$$\|v_a-v_b\|^2 = (\frac{m}{a}-\frac{m}{[a,b]})a+(\frac{m}{b}-\frac{m}{[a,b]})b+\frac{m}{[a,b]}(\sqrt{a}-\sqrt{b})^2$$

$$=2m(1 - \frac{\sqrt{ab}}{[a,b]}) = 2m(1-\frac{(a,b)}{\sqrt{ab}}).$$

This means that after normalization by $2m$ we get the desired embedding. For an embedding in $\mathbb{R}^n$ take the induced embedding in the subspace $span(v_1,...v_n)$.

----------

Another nice embedding straight to a Hilbert space follows from the identity for any natural numbers $a,b$

$$\int_0^1 \psi(at)\psi(bt) dt = \frac{1}{12} \frac{(a,b)^2}{ab}.$$

Where $\psi(t) = x - \lfloor x \rfloor  - \frac{1}{2}$ is the sawtooth function. Therefore for any natural numbers $a,b$

$$\|\psi(at) - \psi(bt)\|_{L^2}^2 = \frac{1}{6}(1-\frac{(a,b)^2}{ab}).$$

So the embedding $\mathbb{N} \hookrightarrow L^2([0,1])$ taking $n \mapsto \psi(nt)$ (also normalizing by $\frac{1}{6}$) preserves this metric!

From the point of view of Fourier series this construction is similar to the previous one, noticing that $\psi(nt)$ only has non-zero Fourier coefficients at the entries divisible by $n$.