Ornstein-Zernike relation from the grand canonical distribution function

Defining the local activity by $z({\bf r})=z\exp[-\beta\psi({\bf r})]$ where $\beta=1/k_BT$, and $k_B$ is the Boltzmann constant. Using those definitions the grand canonical partition function can be written as \begin{eqnarray} \Xi=\sum_N^\infty{1\over N!}\int\dots\int \prod_i^Nz({\b f r}_i)\exp(-\beta U_N){\rm d}{\bf r}_1\dots{\rm d}{\bf r}_N. \end{eqnarray} By functionally-differentiating $\Xi$ with respect to $z({\bf r})$, and utilizing the mathematical theorem concerning the functional derivative, \begin{eqnarray} {\delta z({\bf r})\over{\delta z({\bf r'})}}=\delta({\bf r}-{\bf r'}), \end{eqnarray} we get the following equations with respect to the density pair correlation functions. \begin{eqnarray}\rho({\bf r})={\delta\ln\Xi\over{\delta \ln z({\bf r})}}, \end{eqnarray} \begin{eqnarray} \rho^{(2)}({\bf r,r'})={\delta^2\ln\Xi\over{\delta \ln z({\bf r})\delta\ln z({\bf r'})}}. \end{eqnarray} A relation between $\rho({\bf r})$ and $\rho^{(2)}({\bf r,r'})$ can be obtained after some manipulation as, \begin{eqnarray} {\delta\rho({\bf r})\over{\delta \ln z({\bf r'})}}=\rho^{(2)}({\bf r,r'})-\rho({\bf r})\rho({\bf r'})+\delta({\bf r}-{\bf r'})\rho({\bf r}).\label{deltarho} \end{eqnarray} Now, we define the direct correlation function by an inverse relation of Eq. (\ref{deltarho}), \begin{eqnarray} {\delta \ln z({\bf r})\over{\delta\rho({\bf r'})}}={\delta({\bf r}-{\bf r'})\over{\rho({\bf r'})}} \label{deltalnz}-c({\bf r,r'}). \end{eqnarray} Inserting Eqs. (\ref{deltarho}) and (\ref{deltalnz}) into the chain-rule theorem of functional derivatives, \begin{eqnarray} \int{\delta\rho({\bf r})\over{\delta \ln z({\bf r}^{\prime\prime})}}{\delta \ln z({\bf r}^{\prime\prime})\over{\delta\rho({\bf r'})}}{\rm d}{\bf r}^{\prime\prime}=\delta({\bf r}-{\bf r'}), \end{eqnarray} one get the Ornstein-Zernike Equation. Thus the O-Z equation is, in a sense, a differential form of the partition function.