7

I am stuck with a problem related to simulating a Lennard-Jones system. The Lennard Jones potential is $U(r) = 4\epsilon [ \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} ]$. Hence the force will be $F(r)=-\frac{\partial{U(r)}}{\partial{r}}= 48 \epsilon [ \frac{\sigma^{12}}{r^{13}} - \frac{\sigma^6}{r^7}]$. But to do computer simulation I need to know the force componentwise i.e. all the three components $F_x$, $F_y$, $F_z$ are to be known.

Is the $x$ component simply $F_x = 48 \epsilon [ \frac{\sigma^{12}}{\Delta x^{13}} - \frac{\sigma^6}{\Delta x^7}]$, where $\Delta x$ is the distance between the pair of particles concerned? (or simply $\Delta x =x_i - x_j$)

And let's say that we are calculating the force on the ith particle due to the $j$th particle. How can I express the direction of this force in terms of $x_i$ and $x_j$ ?

Irregular User
  • 3,862
kanayamalakar
  • 191

2 Answers2

7

Just wanted to add (although this is a very old post) that the derivative you have written is not correct, there is a 0.5 factor missing.

$$ F(r) = - \frac{\partial U(r)}{\partial r} = 48 \varepsilon \left[ \frac{\sigma^{12}}{r^{13}}-0.5 \frac{\sigma^6}{r^7} \right]$$

ma7642
  • 73
  • 2
    Thanks for the correction. But note that for numerical simulations (assuming IEEE floats as opposed to arbitrary-precision or symbolic), a more accurate and much faster formula is F/r = 48ε * (σ/r)^6 * ( (σ/r)^6 - 0.5 ) / r^2. About performance, this one doesn’t need r only r^2, i.e. no need for square roots. F/r number is the multiplier to apply to the difference between the coordinates. About accuracy, the subtraction of these high degrees is often a subtraction of very large numbers, precision is degraded due to the nature of binary floating-point representation. – Soonts Feb 09 '22 at 12:54
4

Since the potential is spherically symmetric, the associated force will have only a radial component.

In general, if $U=U(r)$, then $\vec F=-\hat r U'(r)$. We can express this in Cartesian coordinates by using the relationships

$$r=\sqrt{x^2+y^2+z^2}$$

and

$$\hat r =\frac{\hat xx+\hat yy+\hat zz}{\sqrt{x^2+y^2+z^2}}$$

Then, we have

$$\vec F=-\frac{\hat xx+\hat yy+\hat zz}{\sqrt{x^2+y^2+z^2}}\left. \left(\frac{dU(r)}{dr}\right)\right|_{r=\sqrt{x^2+y^2+z^2}}$$

Mark Viola
  • 179,405
  • Doing this I got F(r) as a function of even powers of r along with $(x \hat{x}+ y \hat{y} + z\hat{z})$ term .But can't understand how to write $x \hat{x}$ in terms of $x_i$ and $x_j$ – kanayamalakar Apr 14 '16 at 16:43
  • $\hat x(x_i-x_j)+\hat y(y_i-y_j)+\hat z(z_i-z_j)$ – Mark Viola Apr 14 '16 at 16:46