Question during radius of gyration calculation from XYZ trajectoris of atoms
You don't say exactly what the difference is between your calculation and the LAMMPS calculation, but I'm going to make a guess that you have neglected to account for the periodic boundary conditions properly. The LAMMPS documentation for the compute gyration command says


  The coordinates of an atom contribute to Rg in “unwrapped” form, by using the image flags associated with each atom. See the dump custom command for a discussion of “unwrapped” coordinates. See the Atoms section of the read_data command for a discussion of image flags and how they are set for each atom. You can reset the image flags (e.g. to 0) before invoking this compute by using the set image command.


If you are trying to duplicate this calculation "by hand", you need to "unwrap" the atomic coordinates. This means that, for every molecule, all the atoms in the molecule need to belong to the same "copy" of the molecule, i.e. not to different periodic images of the molecule. You need to ensure this before you even calculate the position of the centre of mass of each molecule.

I'll illustrate one way to do this for a simple linear chain molecule of $N$ atoms. Something similar should work for a more complicated molecule, but you'll need to think about it yourself. I'm going to assume that the atoms are labelled in order $1, 2, 3, \ldots N$, and that $1$ is bonded to $2$, $2$ to $3$, and in general $j$ bonded to $j+1$. Let the atomic positions be $\mathbf{r}_j$. You may want to copy these to a separate working array, to avoid changing the originals.

We'll take the position of atom $1$ as a reference point. In the following loop, start by setting $j=2$.


Compute $\mathbf{d}_{j}=\mathbf{r}_{j}-\mathbf{r}_{j-1}$.
Apply the minimum image correction to $\mathbf{d}_{j}$.
Redefine $\mathbf{r}_{j}=\mathbf{r}_{j-1}+\mathbf{d}_{j}$.
Increment $j\rightarrow j+1$ and return to step 1, unless $j=N$ in which case stop.


Now you should be safe to compute the centre of mass of that molecule, and work out the radius of gyration, in the way you describe in your question. You can repeat the calculation for every molecule independently.

A further point about the LAMMPS compute gyration command is that it not only returns the scalar value $R_g$ (which should correspond to what you describe calculating), it also returns the separate $xx$, $yy$, and $zz$ terms, and also the off-diagonal $xy$, $yz$ and $zx$ terms. This is because the gyration tensor is, like the inertia tensor, a second-rank symmetric tensor. This may be of interest if your polymer configurations are significantly non-spherical in shape. Using these values, you can determine the principal axes of the gyration tensor for each molecule, and the corresponding three principal radii of gyration. It's basically a question of diagonalizing the $3\times3$ symmetric matrix whose components are returned by the LAMMPS command, along with $R_g$. I'm assuming, though, that this is not what you are worried about.
[EDIT: I've just noticed that @porphyrin made this same point in a comment while I was typing in my answer].



EDIT following OP comment.

It seems that "unwrapped" coordinates are already being used for this calculation, in post-processing. So it is possible that my answer does not  explain the difference between the LAMMPS calculation and the one implemented by the OP.

However, different definitions of "unwrapped" may be used for different purposes. In calculating atomic mean-squared displacements as a function of time, to get the diffusion coefficient, it is sufficient to make sure that periodic boundary corrections are never applied to the atomic positions as they evolve in time. However, for that purpose, there is no need to consider which molecule an atom belongs to. In order to compute the molecular centre-of-mass, and from that the radius of gyration $R_g$, the atoms must be grouped together by molecule. This means that, at some stage, a calculation similar to the one described above must be carried out. Maybe this has been done, in producing the "unwrapped" trajectories, maybe not.

So, it is still worth either carrying out the correction described above, or at least checking that bonds between atoms have reasonable lengths in the "unwrapped" trajectories. For example, for a linear chain of $N$ atoms, compute $d_j^2=|\mathbf{r}_j-\mathbf{r}_{j-1}|^2$ for $j=2,3,\ldots N$, and just store the maximum $d_j^2$ in the molecule. Extend this calculation to cover every polymer molecule, and print out the maximum value of $d_j^2$ discovered for the whole set of molecules. How does this compare with the square of the bond length between adjacent atoms, defined by the molecular force field? If there has been no error with the periodic boundaries, the two numbers should at least be similar.