Deriving the 1st law of thermodynamics from statistical mechanics 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?
 A: I think you have to appeal to some sort of physical intuition to show $\frac{\partial E_j}{\partial V_j} = - P.$ In the microcanonical ensemble, we have immediately from your expression for $\rho$,
$$ P = -\beta^{-1} \frac{\partial}{\partial V} \ln Z.$$
Now let us assume $P$ is finite. We pass to the canonical ensemble where $E$ is now a function of $V$ in general. We find
$$ -\frac{\partial}{\partial V} \ln Z = -\frac{1}{Z}\sum_{j} \rho^{-1}_j \sum_i \rho_i \left(-\beta \frac{\partial E_i}{\partial V_j} - \beta P \frac{\partial V_i}{\partial V_j} \right) = \beta Z^{-1}\sum_j\left( -\frac{\partial E_j}{\partial V_j} - P\right).$$
This sum diverges in the thermodynamic limit, unless of course $ \frac{\partial E_j}{\partial V_j} \to -P$ in the right way.
Of course if one assumes Newton's laws (+ other physical assumptions about pressure, etc.) one can also derive $\frac{\partial E}{\partial V} = -P$ from dynamics and then appeal to the microcanonical ensemble. Probably then one could also qualify the convergence $ \frac{\partial E_j}{\partial V_j} \to -P$ given kinetic theory.
EDIT: Using the equivalence of ensembles there are of course many ways to see $\frac{\partial E_j}{\partial V_j} = - P,$ for example by interpreting $P$ as a Lagrange multiplier or by doing a saddle-point approximation on $\rho$.
