# Question during radius of gyration calculation from XYZ trajectoris of atoms

I have molecular dynamic simulation result of 100 polymer molecules (liquid), and I'm trying to calculate the radius of gyration using atomistic mass, molecular mass, and xyz coordinate of atoms. I'm trying to follow what LAMMPS is doing to calculate the radius of gyration $$R_\text{g}$$:

$$R_\text{g}^2 = \frac{1}{M} \sum_i m_i (r_i - r_\text{cm})^2$$

Here is my step:

1) Calculate distance from each atom to COM of molecule, for X, Y, and Z direction.

2) Square the distances (for each X, Y, and Z)

3) Multiply distances (for each X, Y, and Z) by atomic mass for each atom

4) Repeat 1~3 for all atoms of corresponding molecule and sum up the values.

5) Divide the value from 4) by molecular mass

6) Now we have the square of the radius of gyration, but it has 3 components (X,Y,Z). Sum the three components together (they are already squared) and then square root. This would be a Rg of single molecule.

7) Perform 1~6 for all other molecules. (and average in the future)

My calculation result is different from what I get from LAMMPS calculation result. Which part of my script is wrong? I'm not sure what is wrong with my calculation process. Do I need to exclude mass, only consider Ri-Rcom?

• Your method seems to be ok but you could have a typo somewhere in your code .The radius of gyration $R$ is defined for each axis just as is the moment of inertia $I$ and this is defined for each axis $\alpha$ as $I_\alpha=MR_\alpha^2$ where $M$ is the total mass. You might be better calculating the moment of inertia first. If you want a single $R$ then you have to decide what you need, average, average squared etc. – porphyrin Dec 21 '18 at 9:22
• Thinking again, your equations should produce $R_{\alpha,\beta}$ where $\alpha,\beta$ are combinations of x, y and z, i.e (x,x) (x,y), (x,z) etc making 9 in all. These form a square matrix that should be diagonalised to find eigenvalues which are the $R$'s you seek. As mentioned above it is normal to find the moments of inertia first in the matrix method just outlined. – porphyrin Dec 21 '18 at 17:09
• @porphyrin Thank you, Right now, my script calculate dx (Xi-Xcom), dy (Yi-Ycom), and dz (Zi-Zcom), so it seems I only consider Rg-xx, Rg-yy, and Rg-zz. Should I need to consider Rg-xy, Rg-xz, and Rg-yz as well, to calculate overall Rg? – exsonic01 Dec 23 '18 at 2:57
• yes, see my answer below – porphyrin Dec 24 '18 at 11:06

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$$.

1. Compute $$\mathbf{d}_{j}=\mathbf{r}_{j}-\mathbf{r}_{j-1}$$.
2. Apply the minimum image correction to $$\mathbf{d}_{j}$$.
3. Redefine $$\mathbf{r}_{j}=\mathbf{r}_{j-1}+\mathbf{d}_{j}$$.
4. 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.

• thank you so much. My script does not consider off-diagonal parts, only consider xx, yy, and zz. Do I need to include them as well to calculate overall scalar Rg? Regarding LAMMPS, I print out unwrapped coordinates, and all calculations read unwrapped xyz. Problem of LAMMPS is the limit of number of groups, so I can't print out Rg or other data for all molecules, so I need to post process – exsonic01 Dec 23 '18 at 2:59
• No, you don't need to consider the off-diagonal elements in order to compute the scalar $R_g$. Your existing procedure looks OK. If you are already using unwrapped coordinates, I am not sure that my answer deserves to be accepted! It may be that someone else can think of the right answer. – user64968 Dec 23 '18 at 5:21
• I have had a further thought about this, and have edited my answer to suggest a further check that you can carry out if you wish. – user64968 Dec 23 '18 at 7:13

It is easier to start from the beginning and explain how to calculate the moments of inertia. The equations below are a bit different from the one you quote.

If all the atoms are rigidly connected together, the $$k^{th}$$ atom and its velocity vector $$\hat v_k$$, are related to the angular velocity of the molecule $$\hat \omega$$ about the centre of mass as

$$\hat v_k =\hat\omega \times \hat r_k$$

where $$\hat r_k$$ is the position vector from the centre of mass and $$\hat \omega$$ is a vector but does not carry an index. This is because in a rigid body, all atoms move with the same angular velocity. The angular momentum for the $$k^{th}$$ atom is defined as the vector cross product

$$\hat J_k =\hat r_k \times \hat p_k$$

where $$\hat p = m\hat v$$ is the momentum vector and the total angular momentum is the sum over all $$n$$ atoms and is

$$\hat J = \sum_{k=1}^n m_k(\hat r_k \times \hat v_k )=\sum_{k=1}^n m_k(\hat r_k \times (\hat\omega \times \hat v_k )$$

Expanding the triple product into dot products gives

$$\hat J = \sum_{k=1}^n m_k\left( (\hat r_k \cdot \hat r_k)\hat\omega -(\hat r_k \cdot \hat\omega)\hat r_k \right) = \sum_{k=1}^n m_k\left(r_k^2\hat\omega-(\hat r_k \cdot \hat\omega)\hat r_k \right)$$

The vector $$\hat J$$ has components $$x, \,y$$, and $$z$$ so it represents three equations. This can be written as a matrix equation but note that $$r_k^2$$ is a number; it is the perpendicular distance of atom $$k$$ from an axis, $$x,\, y$$ or $$z$$, but $$\hat r_k$$ is the vector $$\hat r_k = (x_k\, y_k \,z_k)$$ describing the position of atom $$k$$.

$$\hat J_{(x,y,z),k}=m_kr_k^2 \begin{bmatrix}\omega_x \\ \omega_y\\ \omega_z \end{bmatrix}-m_k \left( [x_k \,y_k \, z_k] \begin{bmatrix}\omega_x \\ \omega_y\\ \omega_z \end{bmatrix} \right) \begin{bmatrix}x_k\\ y_k\\ z_k \end{bmatrix}$$

The $$x$$ component for the $$k^{th}$$ atom is found by expanding the dot product as $$x_k\omega_x + y_k\omega_y + z_k\omega_z$$ and then multiplying by $$x_km_k$$ and rearranging a little

$$J_{x,k} = m_k\left((r_k^2 - x_k^2)\omega_x - x_ky_k\omega_y - x_kz_k\omega_z \right) \tag{1}$$

There are similar equations for the $$y$$ and $$z$$ direction components.

By comparing coefficients of $$\omega_x$$, and those of $$\omega_y$$ and $$\omega_z$$ the equations the diagonal terms in this matrix are

$$I_{x,x} = \sum_k m_k(r_k^2-x_k^2)$$

and similar equations for $$y$$ and $$z$$.

Because $$r^2=x^2+y^2+z^2$$ we can rewrite $$I_{x,x} =\sum_k m_k (y_k^2 + z_k^2)$$ and similarly for the two other diagonal terms. These terms are called moments of inertia coefficients and cannot be negative as they are the sum of squared terms.

The cross terms $$I_{x,y}$$, for example, are called products of inertia

$$I_{x,y}=-\sum_k m_kx_ky_k, \qquad I_{x,z}=-\sum_k m_kx_kz_k, \qquad I_{y,z}=-\sum_k m_ky_kz_k$$

Equation 1 can be rewritten for each atom $$k$$ using the inertial coefficients

\begin{align} J_x &= I_{xx}\omega_x + I_{xy}\omega_y + I_{xz}\omega_z \\ J_y &= I_{yx}\omega_x + I_{yy}\omega_y + I_{yz}\omega_z\\ J_z &= I_{zx}\omega_x + I_{zy}\omega_y + I_{zz}\omega_z\end{align}

and, of the nine coefficients, only six are different because of symmetry; $$I_{xy} = I_{yx}$$ and so forth. In matrix form these equations are

$$\hat J=\hat I\hat \omega =\begin{bmatrix} I_{xx} & I_{xy} & I_{xz} \\ I_{xy} & I_{yy} & I_{yz} \\ I_{xz} & I_{yz} & I_{zz} \end{bmatrix}\begin{bmatrix} \omega_x \\ \omega_y \\ \omega_z\end{bmatrix}$$

The matrix $$\hat I$$, is also sometimes either called the moment of inertia dyadic or the inertia tensor, but, more importantly, it is symmetrical and Hermitian so has real eigenvalues and orthogonal eigenvectors.

The next step in the calculation is to perform a principal axis transform, which we can view as a rotation of the inertia matrix to remove all the off-diagonal terms that become zero on forming a diagonal matrix, i.e find the eigenvalues of the matrix which produce moments of inertia about the principal axes.

Some python code is given below: it could be refined but should work

# coordinates in xyz[i,j], i is atom, j=0,1,2 = x,y,z ;
# masses in mass(atm[i])
# tmass is total mass of molecule
# coordinates in angstrom, mass in amu

import numpy as np
from numpy import linalg as LA

amu  = 1.6605e-27
angst= 1e-10

com=    [0.0 for i in range(3)]
ss =    np.zeros((n,3),dtype=float )   # n = number atoms,  3 -> x,y,z
for i in range(n):
ss[i,0]=  xyz[i,0]*mass[atm[i]]
ss[i,1]=  xyz[i,1]*mass[atm[i]]
ss[i,2]=  xyz[i,2]*mass[atm[i]]

com = np.sum(ss,0)/tmass
print('{:s} {:10.5g} {:10.5g} {:10.5g}'.format( 'cente of mass', com[0],com[1],com[2]) )

ss = xyz - com
print('CoM based coordinates \n',ss )
r = [0.0 for i in range(n)]
for i in range(n):
r[i]=np.sqrt((xyz[i][0]- com[0])**2+(xyz[i][1]-com[1])**2+(xyz[i][2]-com[2])**2)

Ixx= np.sum( mass[atm[i]] *(r[i]**2  - (ss[i][0])**2 ) for i in range(n)  )
Iyy= np.sum( mass[atm[i]] *(r[i]**2  - (ss[i][1])**2 ) for i in range(n)  )
Izz= np.sum( mass[atm[i]] *(r[i]**2  - (ss[i][2])**2 ) for i in range(n)  )

Ixy= np.sum( -mass[atm[i]] * ss[i][0]*ss[i][1]  for i in range(n)  )
Ixz= np.sum( -mass[atm[i]] * ss[i][0]*ss[i][2]  for i in range(n)  )
Iyz= np.sum( -mass[atm[i]] * ss[i][1]*ss[i][2]  for i in range(n)  )

M= [[Ixx,Ixy,Ixz],[Ixy,Iyy,Iyz],[Ixz,Iyz,Izz]]
eigvals = LA.eigh(M)                             # eigenvalues

print('\n eigenvalues\n', eigvals[0])
print('\n moments of inertia kg.m^2', eigvals[0]*amu*angst**2)
print('radius of gyration /m', np.sqrt(eigvals[0]*amu*angst**2/(tmass*amu) ) )