Two-body correlation function computation

I'm studying how to compute excess entropy in molecular dynamics (MD). I've found it is needed to compute the two-body correlation function (neglecting high-order terms), the details can be found, for example, in this article.

So the definition of correlation function (CF for short) is
##C(t, r,t',r')=\langle X(t,r)Y(t',r')\rangle##
where angle brackets mean averaging.

First question: is the averaging performed by time or ensemble (by all the atoms in the system)?

Second: for computing the CF, do I need to compute it in the stationary process? I mean, do I need to simulate a steady-state system in MD (probably to perform time-averaging) or can CF be found from one time point (using atom coordinates and velocities in a specific time moment)?

Third: If, for example, I want to compute CF for relative distance ##\vec r=\vec r_2 - \vec r_1##, where ##\vec r_1, \vec r_2## are the absolute positions of two atoms, what will be ##X## and ##Y##?I'm sorry if I've written something unclear, I'm always ready to clarify the question, and I'd be happy for any help.

P.S. My goal is to calculate the free Gibbs energy of bunch of atoms via computing entropy or configuration integral (partition function).

It seems like I have understood the correlation function idea.

The averaging type depends on property we have to study, so for relative distance averaging is performed by all particles. And in case of correlation for relative distance CF is turned into radial distribution function, so to calculate it we have to take one particle, calculate distances to other particles, make a histogram with given $\Delta r$, repeat these steps for all particles, calculate the average bin heights and normalize obtained distribution. The same is for angular position CF $$g(\phi_{pos}, \theta_{pos})$$.