Here's an alternative proof. It seems more elementary to me, but that might just be because I don't see the connections to the results used in David's proof and am reinventing the wheel. If there are such connections, I'd appreciate if someone could point them out in the comments. My proof will use David's result that $|\sum a_k \zeta^{k}|=1$ implies $|\sum a_k \zeta^{kr}|=1$ for $r$ relatively prime to $n$; the only "sophisticated" result he used for that was the irreducibility of the cyclotomic polynomials (and implicitly unique factorization).
First, note that the question has changed and the answer to the new question is "no". As I wrote in the comments, two cube roots of unity sum to a square root of unity, so the restriction to $n$-th roots of unity makes the statement false.
The true statement in the original question was "If the sum of some $n$th roots of unity has magnitude $1$, this implies that this sum is some root of unity as well" (where the same $n$-th root of unity may occur multiple times in the sum, to make this equivalent to the integer linear combination formulation). This is equivalent to the same statement with "$n$-th" deleted, since for every set of roots of unity there is a common denominator $n$ for which they are all $n$-th roots.
Now for the proof, consider the free vector space on the set of the $n$-th roots of unity. Any linear combination of $n$-th roots of unity acts on this vector space by multiplication. The corresponding matrix with respect to the canonical basis is a circulant matrix with the first column given by the coefficients $a_i$ of the linear combination. The eigenvectors of this circulant matrix are the Fourier modes $v_r=(1,\omega^r,\omega^{2r},\dotsc$), with $\omega=\mathrm e^{2\pi\mathrm i/n}$, and the corresponding eigenvalues are the Fourier transforms of the coefficients, $\sum a_k \omega^{-kr}$.
Now consider the natural map $f$ from this vector space to $\mathbb C$ sending each vector to the corresponding linear combination of the roots of unity. All but one of the Fourier modes are in the kernel of this map. Only the $r=-1$ mode isn't, and the corresponding eigenvalue is our number $\alpha=\sum a_k \omega^k$.
Now let's start with the vector $(1,0,\dotsc,0)$, corresponding to unity, and successively multiply by the circulant matrix. This initial vector has equal components $1/n$ of each Fourier mode. Consider what happens if we add the components corresponding to the primitive roots for some divisor $d\mid n$ together. If $d\neq n$, all of these are in the kernel of $f$. The entries in their sum are Ramanujan's sums, which are integers. Thus, by dropping these components, we don't change the complex number that the vector maps to, and we don't change the fact that all entries in the vector are integer multiples of $1/n$.
Now all that's left is the sum of the components corresponding to the primitive $n$-th roots. We know that the eigenvalue for $r=-1$ has magnitude $1$, and by David's result that implies that all the others also have magnitude $1$. So we have a vector whose entries are integer multiples of $1/n$ and we keep multiplying it by an integer matrix whose eigenvalues (as far as this vector is concerned) have magnitude $1$. It follows that there are only finitely many values this vector can take, so it has to eventually get back to where it started, which implies that $\alpha$ is a root of unity.