# How do I calculate the radial distribution function from the centre of mass of unit cell in GROMACS?

I am simulating a droplet of water with GROMACS. The droplet consists of about 200 TIP3P water molecules and is placed near the center of a large, mostly empty box (so that it behaves like a free droplet in a vacuum, in spite of periodic boundary conditions. I had to do this as GROMACS2021 does not allow a simulation without PBC).

Now, I have to calculate the radial distribution function of the oxygen atoms and hydrogen atoms from the center of mass of the whole unit cell. However, I am not sure how to do this with gmx rdf. The documentation on the RDF module is minimal, and not clear at all.

I tried to run this:

gmx rdf -f md.xtc -s md.tpr -n index.ndx -o com-O.xvg -ref "com of group 0" -sel 3


In the index.ndx file, group 0 is System, group 3 is O atoms. However, this gave a graph that looks like this:

The graph is very squiggly, and I expected a smoother graph.

I have also tried to use -ref 0 -selrpos dyn_mol_com but that gives a graph even more uneven.

So my question is what commands should I use with gmx rdf to get the rdf, and whether my graph looks ok or not?

Note: I have already gone through these two pages: gmx rdf and selection syntax and usage, but I cannot understand what they mean. Please don't redirect me to those pages!

• +1. I don't have the answer for you, but I thought I'd link some related questions here. About RDFs in GROMACS: mattermodeling.stackexchange.com/q/3976/5, and about RDFs in general: mattermodeling.stackexchange.com/q/3852/5, mattermodeling.stackexchange.com/q/3624/5, mattermodeling.stackexchange.com/q/90/5. – Nike Dattani Mar 15 at 23:01
• Are you open to using other packages, like mdtraj or mdanalysis (both python)? – lewiso1 Mar 21 at 23:09
• @lewiso1 Do you mean other packages just for trajectory analysis, or the simulation itself? I can use any software for trajectory analysis, but the simulation has to be ran in gromacs. – Shoubhik R Maiti Mar 22 at 11:24
• I added an answer using numpy, but just realised you can probably get rid of the squiggliness by just reducing the binwidth. The droplet only has a radius of 1nm, which is quite small, so currently the RDF is way higher resolution than what you need and there's a lot of noise. Try adding -bin 0.15 to your gmx rdf – lewiso1 Mar 23 at 0:26

Here's an approach calculating it by hand with numpy. The coordinates are loaded from the trajectory using mdtraj in python:

#imports:
import mdtraj as md
import numpy as np
import matplotlib.pyplot as plt

traj = md.load('traj.dcd', top='solvated.pdb')
##You would write this:
##traj = md.load_xtc('md.xtc', top='top.gro')

#select oxygen atom indices:
oxy = traj.top.select('name O')

#calculate masses of all atoms:
masses = np.array([i.element.mass for i in traj.topology.atoms])

#iterate through all frames and save the distances of each atom to the centre of mass:
d = list()
for frame in traj.xyz:
oxy_coords = frame[oxy]
#COM of whole system:
centre_of_mass = np.average(frame, axis=0, weights=masses)
#COM of just the oxygens:
#centre_of_mass = oxy_coords.mean(0)
vectors = (oxy_coords - centre_of_mass)
distances = np.linalg.norm(vectors, axis=1)
d += list(distances)

##Take a histogram:
counts, lengths = np.histogram(d,bins=10)
#normalize 'counts' by averaging across the number of trajectory frames:
counts = counts / len(traj)
#calculate the volume of each spherical shell:
shell_volumes = 4/3*np.pi* ( lengths[1:]**3 - lengths[:-1]**3 )
#normalize 'counts' by the volume of each shell:
counts = counts / shell_volumes

#plot:
plt.plot( (lengths[:-1]+lengths[1:])/2,
counts , '-o' )
plt.ylabel('Number density per nm^3')
plt.xlabel('r (nm)')
plt.axhline(33.6, label='Waters per cubic nm\nat room temp and pressure', c='k')
plt.legend()


I tried this on an example trajectory and it returns the graph below, which has number density ~33, which is what you would expect at room temperature and pressure. Note that I extracted a water droplet from a simulation of a periodic liquid system with a barostat, so your result will be different if you used a water droplet as the simulation system.

• Thanks! One thing I noticed is the line centre_of_mass = oxy_coords.mean(0). I assume this means the centre of mass is being estimated by the oxygen atoms only? How do I also consider the hydrogen atoms because I want to use the precise centre of mass of the unit cell? – Shoubhik R Maiti Mar 23 at 9:59
• Sure - I edited the code to do this. The difference is you take an array of masses, and then when calculating centre_of_mass use np.average on all coordinates to find the mean coordinate, weighted by masses – lewiso1 Mar 25 at 2:17