Here is a method.
Let $\mathrm{V}$ be the bilinear space $(\mathbf Q^n, b)$ with Gram Matrix $G$ in the canonical basis.
Let $M_0$ be the lattice $\mathbf Z^n$ on $\mathbf Q^n$. Find a primitive vector $v_0\in M_0$ such that $b(v_0,v_0)$ is a square, say $b(v_0,v_0)=m^2$ in $M_0$ (note that this can be done quite efficiently : use LLL to get a basis of vectors with small norms, take the two* first vectors, study the lattice they generate).
( * Edit 1: in fact you need to take the first 4 vectors ... indeed, there exists (many) $3\times 3$ positive definite Gram matrices that don't represent any square over $\mathbf Q$ ... The simplest example is $7.I_3$.)
( Edit 2: If you don't use this trick, you can arrange so as to find a $v_0$ such that $m=2^s$ ... this will allow you to find an orthonormal basis inside $\mathbf Z[\frac 12]^n$, whose existence is predicted by the Strong Approximation Theorem.)
Let $K_1$ be the lattice made of vectors $w\in M_0$ such that $(v,w)\equiv 0$ mod. $m$.
Let $L_1$ be the lattice $K_1+\mathbf Z\cdot\frac{v_0}{m}$. Then $L_1$ is integral, unimodular, and has a vector of length $1$, namely $w_0:=\frac{v_0}{m}$.
Let $M_1$ be the orthogonal complement of $w_0$ in $L_1$. Then $M_1$ is an $n-1$-dimensional unimodular lattice (contained in the $n-1$-dimensional subspace $<w_0>^\perp$ of $\mathrm{V}$).
Go on to find $v_1\in M_1$, construct $M_2$ and so on ... you're done (the matrix $M$ you were looking for is the inverse of the matrix whose columns are $(w_0,\dots,w_{n-1})$.
[The classification of unimodular lattices shows that at least when $n\leq 7$, you can always arrange so as to have $(v_0,v_0)=1$ (so that $M_1\subset M_0$), and for $n\leq 24$, you can find some $v_0$ with $(v_0,v_0)=4$.]