I'm attempting to derive the general 1st law of thermodynamics
$$\mathrm{d}\langle E\rangle = T\,\mathrm{d}S + \sum_{i = 1}^\ell P_i\,\mathrm{d}\langle X_i\rangle,$$
where $\langle E\rangle$ is the average energy of the system, $S$ is the entropy of the system, $\langle X_i\rangle$ is the average value of the $i^\text{th}$ extensive variable of the system (say volume or particle number), and $P_i$ is the $i^\text{th}$ intensive variable of the system (say pressure or chemical potential). In particular, suppose that we are considering a system which may exchange only energy and volume with its environment (the general case in which the system can exchange arbitrary extensive parameters with its environment should generalize easily). Starting from Liouville's theorem, it can be shown that the phase space distribution function $\rho$ describing a macroscopic system must take the form
$$\rho(E,V) = \frac{1}{Z}\exp(-\beta E - \beta P V),$$
where $\beta$, $P\beta$, and $\mu\beta$, are constants (which we will later interpret as the inverse temperature, pressure, and chemical potential), and where $\frac{1}{Z}$ is a normalization constant. Define the entropy as
$$S = -k_B \sum_j \rho_j \ln(\rho_j),$$
where $\rho_j = \rho(E_j,V_j)$. Taking the differential gives
$$\mathrm{d}S = -k_B\sum_j\left(\ln(\rho_j)+1\right)\mathrm{d}\rho_j = -k_B\sum_j\left(-\ln(Z)-\beta E_j - \beta P V_j+1\right)\mathrm{d}\rho_j.$$
The quantity $\sum_j\mathrm{d}\rho_j$ vanishes due to normalization, and so we have that
$$\mathrm{d}S = k_B\beta\sum_j\left(E_j + P V_j\right)\mathrm{d}\rho_j.$$
Rearranging, we obtain
$$0 = \frac{1}{k_B\beta}\mathrm{d}S - \sum_j\left(E_j + P V_j\right)\mathrm{d}\rho_j.$$
We now add $\mathrm{d}\langle E\rangle = \sum_jE_j\mathrm{d}\rho_j + \sum_j \rho_j\mathrm{d}E_j$ and $P\mathrm{d}\langle V\rangle = P\sum_jV_j\mathrm{d}\rho_j + P\sum_j \rho_j\mathrm{d}V_j$ to both sides of this equation, giving
$$\mathrm{d}\langle E\rangle + P\mathrm{d}\langle V\rangle = \frac{1}{k_B\beta}\mathrm{d}S + \sum_j \rho_j\mathrm{d}E_j + P\sum_j \rho_j\mathrm{d}V_j.$$
Using $\mathrm{d}E_j = \frac{\partial E_j}{\partial V_j}\mathrm{d}V_j$, we find that
$$\mathrm{d}\langle E\rangle = \frac{1}{k_B\beta}\mathrm{d}S - P\mathrm{d}\langle V\rangle + \sum_j \rho_j\left(\frac{\partial E_j}{\partial V_j} + P\right)\mathrm{d}V_j.$$
Unfortunately, I am stuck here. This is very close to the first law, but I simply cannot figure out what to do with the third term on the RHS. I'm tempted to say that $\frac{\partial E_j}{\partial V_j} = -P$ and be done with it, but I'm not sure that that is valid since pressure is a macroscopic quantity, and not simply equal to $\frac{\partial E_j}{\partial V_j}$ in the $j^\text{th}$ microstate. Any advice on how to proceed?