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

Question

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!

Answer 1

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

#load trajectory:
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.