Local averages on arbitrary three-dimensional grid

GROMACS version: 2024.1
GROMACS modification: No

My aim is to compute local measurements of the compressibility of water around a biomolecule (e.g. a small protein to start with). A way to do it is to define a spatial three-dimensional grid and then compute histograms of suitable quantities linked to the local number of particle. In particular, I am interested in working in the isobaric ensemble (so pressure coupling, e.g. MTTK) and calculate the following correlator

C(x,y,z) = ‹ rho(x,y,z) V› -‹rho(x,y,z)› ‹V›

where rho(x,yz) is the instantaneous density of water (e.g. Oxygen atoms) at r=(x,y,z) and V is the instantaneous volume of the global simulation box, which in the isobaric ensemble fluctuates. ‹› indicate a time average over the equilibrium steady state. Typically, one needs many tens of thousands of independent profiles to produce reasonably smooth local profiles (this from tests in LAMMPS).

My questions are :

  1. is there a recommended way to implement this local averaging on an arbitrary grid as an additional feature in GROMACS ?
  2. if modifying GROMACS source code is complex, is there a recommended workflow that you would suggest that would avoid producing too large files?