Derivation of Pressure Correction Term in Lennard-Jones Fluid Simulation

Question

Context

Hey everyone, I'm working on a molecular dynamics assignment where I'm simulating a Lennard-Jones fluid in the Gibbs ensemble. The potential is given by equation (1): $$ U(r) = 4\epsilon \left( \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right) \tag{1} $$ where we use a cut-off distance $r_c$ to save on computational resources, essentially ignoring interactions that occur beyond a certain point.

Because of this, there then needs to be a correction term to find the total pressure accurately. The equation for this correction is given as equation (2): $$ p_{\text{corr}} = -\frac{1}{6} \left(\frac{N}{V}\right)^{2} \int_{r_{c}}^{\infty} r \frac{d}{d r}\left(-\frac{4 \epsilon \sigma^{6}}{r^{6}}\right) 4 \pi r^{2} d r \tag{2} $$

This simplifies to equation (3): $$ -\frac{16 \pi N^{2} \epsilon \sigma^{6}}{3 V^{2} r_{c}^{3}} \tag{3} $$

Confusion

I get that we're approximating the potential to just its attractive part and then taking its derivative to find the force, as well as assuming the pair distribution function $g(r) = 1$ past the cut-off. But here's where I'm stuck: I can't quite grasp where this equation comes from and how it sums up the attractive force (at long range) for every pair in the system.

As far as what I've "tried", I've been looking for resources that could help me understand. These are the relevant concepts I've seen so far:

Virial Expansion

There is this equation of state given by equation (4): $$ P = \rho k_B T \left(1 + B_2 \rho + B_3 \rho^2 + \ldots \right) \tag{4} $$ where the second virial coefficient is given by equation (5): $$ B_{2}(T)= - \dfrac{1}{2} \int \left( \exp\left(-\dfrac{\Phi_{12}({\mathbf r})}{k_BT}\right) -1 \right) 4 \pi r^2 dr \tag{5} $$

I have a hunch that the $ \rho^2 $-term in equation (4) might relate or equate to $ p_{corr} $, which would explain how we get the squared number density in equation (2). But I can't quite reconcile it with equation (2), so I'm lacking a clear connection between the two. To get them to match I imagine we would have to ignore the linear term and orders above 2, which I don't quite understand the justification for. I'm also not sure how the $ k_B T $ would be involved since it's not in the integrand of equation (2). Is the lack of temperature dependence in (2) due to some sort of approximation? I'm not sure.

Request for Help

I'm juggling courses in statistical mechanics and computational physics at the moment. We haven't touched on this in my stat mech class yet, and I'm really keen to figure out how to derive $ p_{corr} $ and understand its construction. Would really appreciate any insights on:

• How the power series in equation (4) might relate to my special case (2)
• If it does, how I might derive $ p_{corr} $ from it
• If it doesn't, what other resources do I look up
• How from (2) we indeed get the (attractive) interaction force from every pair in the system.

(Just to note, this isn't for a homework problem or anything; I'm just looking for some guidance to deepen my understanding on an equation that was given, so I apologize in advance if this isn't the type of question that should be posted.)

Answer 1

Using the $g(r)$ you can express any average quantity (that depends on the particle-particle distance $r_{ij}$) in integral form:

$$ \left \langle \sum_{i \neq j} a(r_{ij}) \right \rangle = \frac{N^2}{2V} \int_0^{\infty} a(r) g(r) 4 \pi r^2 dr $$

since the sum has $N(N-1) / 2 \approx N^2 / 2$ terms.

Recalling that the pressure can be written as

$$ P = \rho k_B T + \frac{1}{3V} \left\langle \sum_{i \neq j} \vec{f}_{ij} \cdot \vec{r}_{ij} \right\rangle, $$

for central forces $\vec{f} \cdot \vec{r} = rf$ and therefore one finds

$$ \frac{1}{3V} \left\langle \sum_{i \neq j} \vec{f}_{ij} \cdot \vec{r}_{ij} \right\rangle = \frac{N^2}{6V^2} \int_0^\infty r f(r) g(r) 4 \pi r^2 dr. $$

The difference between the "real" pressure and the one we compute in simulations (for which $f(r) = 0$ for $r > r_c$) is then

$$ P_{\rm real} - P_{\rm cut} = \frac{N^2}{6V^2} \int_0^\infty r f(r) g(r) 4 \pi r^2 dr - \frac{N^2}{6V^2} \int_0^{r_c} r f(r) g(r) 4 \pi r^2 dr = \frac{N^2}{6V^2} \int_{r_c}^\infty r f(r) g(r) 4 \pi r^2 dr $$