Questions regarding a derivation in Understanding Molecular Simulation: From Algorithms to Applications, 3rd ed, pp 141-143

Question

Recently, I've been learning molecular dynamics using the third edition of the book, Understanding molecular simulation: from algorithms to applications. During reading your book, I have encountered some confusion from page 141 to page 143 in section 5.1.7.2, and your help shall be very valuable.

Firstly, on page 142 between formula (5.1.47) and (5.1.48), the authors replaced $\hat{\mathbf{r}} \cdot \nabla_{\mathbf{r}}$ with $-\hat{\mathbf{r}_j} \cdot \nabla_{\mathbf{r}_j}$, using an assumption that $\hat{\mathbf{r}} \cdot \nabla_{\mathbf{r}} \delta(\mathbf{r}-\mathbf{r}_j) = -\hat{\mathbf{r}_j} \cdot \nabla_{\mathbf{r}_j} \delta(\mathbf{r}-\mathbf{r}_j)$. This assumption is not true for a generic central symmetric function centered on $\mathbf{r}_j$, so how can I show that it is true for Dirac-$\delta$ functions?

Secondly, in equation (5.1.52) of the book, the authors attempted to regenerate the results of Borgis et al., equation (6). However, in my opinion, the book's result is not evidently the equivalent to the result of Borgis et al. The book wrote in angled bracket $F^{(r)}_j-F^{(r)}_i$, where $F^{(r)}_j$ represents $\hat{\mathbf{r}}\cdot \mathbf{F}_j$, the projection of the force acting on particle $j$ over unit vector $\hat{\mathbf{r}}$. In contrast, the paper of Borgis et al. wrote $(\mathbf{F}_j - \mathbf{F}_i) \cdot \hat{\mathbf{r}_{ij}}$, which is not equal to $F^{(r)}_j-F^{(r)}_i$ as given by the authors (at least not straightforward).

A brief scan of the aforementioned pages are attached. Any suggestion will be highly appreciated!

Reproduction of the relevant derivation from pp141 to 143 in the book (as the original text, be careful as any possible mistake in the original text is not corrected): The value of the radial distribution function at a distance $r$ from a reference particle is equal to the angular average of $\rho(\mathbf{r})/\rho$: $$g(r)=\frac{1}{\rho}\int d\hat{\mathbf{r}} \left< \rho(\mathbf{r})\right>_{N-1}=\frac{1}{\rho} \int d\hat{\mathbf{r}} \left<\sum_{j\neq i}\delta(\mathbf{r}-\mathbf{r}_j) \right>_{N-1}$$ (5.1.44).

Where $N$ is the total number of particles in the system, $\rho$ denotes the average number density and $\mathbf{r}_j$ is the distance of particle $j$ from the origin. \hat{\mathbf{r}} is the unit vector in the direction of $\mathbf{r}$. For simplicity, we have written down the expression for $g(r)$ for a given particle $i$, and hence the sum of $j \neq i$ is keeping $i$ fixed, but in practice the expression should be averaged over all equivalent $i$. The angular bracket denotes the thermal average: $$\left<\cdots\right>_{N-1}=\frac{\int d\mathbf{r}^{N-1} e^{-\beta U(\mathbf{r}^N)(\cdots)}}{\int d\mathbf{r}^{N-1} e^{-\beta U(\mathbf{r}^N)}}$$ (5.1.45).

We can now write $$(\frac{\partial g(r)}{\partial r}) = \frac{1}{\rho} \frac{\partial}{\partial r} \int d\hat{\mathbf{r}} \left<\sum_{j \neq i} \delta(\mathbf{r}-\mathbf{r}_j)\right>$$ (5.1.46)

The only term that depends on $r$ is the $\delta$-function.Thus: $$(\frac{\partial g(r)}{\partial r}) = \frac{1}{\rho} \int d\hat{\mathbf{r}} \left<\sum_{j \neq i} \hat{\mathbf{r}} \cdot \nabla_{\mathbf{r}} \delta(\mathbf{r}-\mathbf{r}_j)\right>$$ (5.1.47)

Since the argument of the $\delta$-function is $\mathbf{r}-\mathbf{r}_j$, we can replace $\hat{\mathbf{r}} \cdot \nabla_{\mathbf{r}}$ by $-\hat{\mathbf{r}}_j \cdot \nabla_{\mathbf{r}_j}$ and perform a partial integration:

\begin{align} (\frac{\partial g(r)}{\partial r}) &= \frac{-1}{\rho} \frac{\int d\hat{\mathbf{r}} \int d\mathbf{r}^{N-1} e^{-\beta U(\mathbf{r}^N)} \sum_{j \neq i} \hat{\mathbf{r}} \cdot \nabla_{\mathbf{r}} \delta(\mathbf{r}-\mathbf{r}_j)}{\int d\mathbf{r}^{N-1} e^{-\beta U(\mathbf{r}^N)}}\\ &=\frac{-\beta}{\rho} \frac{\int d\hat{\mathbf{r}} \int d\mathbf{r}^{N-1} e^{-\beta U(\mathbf{r}^N)} \sum_{j \neq i} \delta(\mathbf{r}-\mathbf{r}_j) \hat{\mathbf{r}_j} \cdot \nabla_{\mathbf{r}_j} U(\mathbf{r}^N)}{\int d\mathbf{r}^{N-1} e^{-\beta U(\mathbf{r}^N)}}\\ &=\frac{\beta}{\rho} \int d\hat{\mathbf{r}} \left<\sum_{j \neq i} \delta(\mathbf{r}-\mathbf{r}_j) \hat{\mathbf{r}}_j \cdot \mathbf{F}_j(\mathbf{r}^N)\right>_{N-1} \end{align} (5.1.48)

where $\hat{\mathbf{r}} \cdot \mathbf{F}_j = F_j^{(r)}$ denotes the force on particle j in the radial direction. We can now integrate with respect to r

\begin{align} g(r) &= g(r=0) + \frac{\beta}{\rho} \int_0^r dr' \int d\hat{\mathbf{r'}} \left<\sum_{j \neq i} \delta(\mathbf{r}-\mathbf{r}_j)F_j^(r)(\mathbf{r}^N)\right>_{N-1}\\ &= g(r=0) + \frac{\beta}{\rho} \int_{r'<r} d\mathbf{r}' \left< \frac{\sum_{j \neq i}\delta(\mathbf{r}-\mathbf{r}_j)F_j^{(r)}(\mathbf{r}^N)}{4\pi r'^2} \right>_{N-1}\\ &= g(r=0) + \frac{\beta}{\rho} \sum_j \left< \theta(r-r_j) \frac{F_j^{(r)}(\mathbf{r}^N)}{4\pi r_j^2} \right>_{N-1} \end{align} (5.1.49)

where $\theta$ denotes the Heaviside step function. To make a connection to the results of Borgis et al., we note that in a homogeneous system, all particles $i$ of the same species are equivalent. We can therefore write:

$$g(r)=g(r=0)+\frac{\beta}{N\rho} \sum_{i=1}^{N} \sum_{j \neq i} \left< \theta(r-r_j) \frac{F_j^{(r)}(\mathbf{r}^N)}{4\pi r_j^2} \right>_{N-1}$$

But $i$ and $j$ are just dummy indices. So we obtain the same expression for $g(r)$ by permuting $i$ and $j$, except that if $\hat{\mathbf{r}}=\hat{\mathbf{r}}_{ij}$, then $\hat{\mathbf{r}}=-\hat{\mathbf{r}}_{ji}$. Adding the two equivalent expressions and dividing by 2, we get

$$g(r)=g(r=0)+\frac{\beta}{2N\rho} \sum_{i=1}^{N} \sum_{j \neq i} \left< \theta(r-r_{ij}) \frac{F_j^{(r)} (\mathbf{r}^N) - F_i^{(r)} (\mathbf{r}^N)}{4\pi r_{ij}^2} \right>_{N-1}$$ (5.1.50)

Eq (5.1.50) is equivalent to the results of Borgis et al.

Answer 1

For the first question, note that if $f$ is a differentiable function then

\begin{equation} \begin{split} -\mathbf y \cdot \nabla_{\mathbf y} f(\mathbf x - \mathbf y) &= -\sum_i y^i \frac{\partial}{\partial y^i} f(\mathbf x - \mathbf y) \\ &= -\sum_{i,j} y^i (D_j f)(\mathbf x - \mathbf y) \frac{\partial}{\partial y^i} (x^j - y^j) \\ &=+\sum_{i,j} y^i (D_j f)(\mathbf x - \mathbf y) \delta_{ij}\\ &= \sum_i y^i (D_i f)(\mathbf x - \mathbf y) \\ &= \mathbf y \cdot \nabla_{\mathbf x} f(\mathbf x - \mathbf y) \end{split} \end{equation}

where $\{x^i\}$ and $\{y^i\}$ are the cartesian components of $\mathbf x$ and $\mathbf y$, respectively, and $D_i f$ denotes the first partial derivative of $f$ with respect to its $i$th argument. Now, suppose that $f$ is somehow zero everywhere except at the origin, i.e. $f(\mathbf x - \mathbf y) = 0$ except when $\mathbf x = \mathbf y$. Certainly then $\nabla f$ is zero away from the origin as well. Thus we can safely exchange $\mathbf x$ and $\mathbf y$ in the above to get

$$-\mathbf y \cdot \nabla_{\mathbf y} f(\mathbf x - \mathbf y) = \mathbf x \cdot \nabla_{\mathbf x} f(\mathbf x - \mathbf y).$$

Finally, the premultipliers $\mathbf x$ and $\mathbf y$ can be replaced by their unit vector counterparts by the same kind of trick.

Of course no differentiable function can be zero everywhere except at the origin (that would be a contradiction), but the delta function behaves as though it was such a function. (It is actually a distribution.)

Answer 2

In case anyone on this forum is still interested, I have written to Prof. Daan Frenkel who is an author of this book, and have received his correction note on the issues mentioned above. Here I give my special thanks to Prof. Daan Frenkel for putting his effort into answering a beginner's questions and into improvising his book! According to Prof Frenkel, there was a few typos and a confusing statement in the original book, but the validity of the final result should not be affected.

Here is his reply:

The following sentence was incorrect: As the argument of the $\delta$-function is $\mathbf{r} - \mathbf{r}_j$, we can replace $\hat{\mathbf{r}} \cdot\nabla_{\mathbf{r}}$ by $-\hat{\mathbf{r}_j} \cdot\nabla_{\mathbf{r}_j}$ and perform a partial integration. This sentence is wrong, because the order is actually the following:

\begin{aligned} & \left(\frac{\partial g(r)}{\partial r}\right)=\frac{1}{\rho} \frac{\int d \hat{\mathbf{r}} \int d \mathbf{r}^{N-1} \mathrm{e}^{-\beta U\left(\mathbf{r}^N\right)} \sum_{j \neq i} \hat{\mathbf{r}} \cdot \nabla_{\mathbf{r}} \delta\left(\mathbf{r}-\mathbf{r}_j\right)}{\int d \mathbf{r}^{N-1} \mathrm{e}^{-\beta U\left(\mathbf{r}^{\mathbb{N}}\right)}} \\ & =\frac{-1}{\rho} \frac{\int \mathrm{d} \hat{\mathbf{r}} \int \mathrm{d} \mathbf{r}^{\mathrm{N}-1} \mathrm{e}^{-\beta U\left(\mathbf{r}^{\mathbb{N}}\right)} \sum_{\mathbf{j} \neq \mathbf{i}} \hat{\mathbf{r}} \cdot \nabla_{\mathbf{r}_j} \delta\left(\mathbf{r}-\mathbf{r}_j\right)}{\int \mathrm{d} \mathbf{r}^{N-1} \mathrm{e}^{-\beta U\left(\mathrm{r}^{\mathbb{N}}\right)}} \\ & =\frac{-\beta}{\rho} \frac{\int \mathrm{d} \hat{\mathbf{r}} \int \mathrm{d} \mathbf{r}^{\mathrm{N}-1} e^{-\beta U\left(\mathbf{r}^{\mathbb{N}}\right)} \sum_{j \neq i} \delta\left(\mathbf{r}-\mathbf{r}_j\right) \hat{\mathbf{r}} \cdot \nabla_{\mathbf{r}_j} \mathrm{U}\left(\mathbf{r}^{\mathrm{N}}\right)}{\int \mathrm{d} \mathbf{r}^{\mathrm{N}-1} e^{-\beta U\left(\mathbf{r}^N\right)}} \\ & =\frac{-\beta}{\rho} \frac{\int \mathrm{d} \hat{\mathbf{r}} \int \mathrm{d} \mathbf{r}^{\mathrm{N}-1} e^{-\beta \mathrm{U}\left(\mathbf{r}^{\mathbb{N}}\right)} \sum_{j \neq i} \delta\left(\mathbf{r}-\mathbf{r}_j\right) \hat{\mathbf{r}}_j \cdot \nabla_{\mathbf{r}_j} \mathrm{U}\left(\mathbf{r}^{\mathrm{N}}\right)}{\int \mathrm{d} \mathbf{r}^{\mathrm{N}-1} e^{-\beta U\left(\mathbf{r}^{\mathbb{N}}\right)}} \\ & =\frac{\beta}{\rho} \int d \hat{\mathbf{r}}\left\langle\sum_{j \neq i} \delta\left(\mathbf{r}-\mathbf{r}_j\right) \hat{\mathbf{r}}_j \cdot \mathbf{F}_j\left(\mathbf{r}^N\right)\right\rangle_{N-1} \\ & \end{aligned}

Also, the sentence "where $\hat{\mathbf{r}} \cdot \mathbf{F}_j = F^{(r)}_j$ denotes the force on particle j in the radial direction" is wrong, because $F^{(r)}_j$ actually represents $\hat{\mathbf{r}_j} \cdot \mathbf{F}_j$ according to our previous derivations. This should also answer your second question, because in the beginning of the derivation, we have centered the reference frame on an arbitrary particle $i$, and therefore the $\hat{\mathbf{r}}_j$ here is actually $\hat{\mathbf{r}}_{ij}$, which points on the direction from particle $i$ to partical $j$. Our result is thus identical to that of Borgis et al.