It should be fairly easy to write a tcl script to do this. Something on the lines of :

1. Set bulk1 as the oxygen of water and bulk2 as the atom corresponding to the center of mass of your osmolyte molecule.

2. Calculate the number of bulk water and bulk solvent molecules using

llength [lsort -unique [$bulk1 get index]]

3. Let r be the distance from the center of mass of the protein.

Likewise, define vicinal1 and vicinal2 as the oxygens and com of osmolyte within $r respectively.

4. Use a nested loop to iterate over varying values of $r and the frames in your trajectory, and calculate the number of vicinal1 and vicinal2 molecules as in step 2.

5. Divide $bulk2 by $bulk1, multiply by $vicinal1 and subtract from $vicinal2 (as per the definition of preferential interaction coefficient)

6. Average over the number of frames to get a time average of the value and plot it as a function of $r.

Thank you for the reply. I have calculated it by using gromacs command gmx select. Is it right to do so? As when i have plotted the graph it continuously goes on increasing rather than showing a dip in the curve.