Determine Forces from a Bond Angle Potential

Question

I am looking for a physics based proof for how one applies forces to each of the three atoms in a bond angle potential. I know the general principle that $F = -\nabla \Phi$ if $\Phi$ is your potential energy function. I am wondering how one goes from such a definition to the force on each atom and how to prove that those forces satisfy conservation of energy and momentum (or any other relevant physical laws).

For some more details, here is how LAMMPs applies forces on atoms from a simple harmonic bond angle potential:

$\Phi(\theta) = k(\theta - \theta_{0})^{2}$

Consider this image:

Then $\Phi = k(\arccos(\frac{\mathbf{r_{1}\cdot r_{2}}}{r_{1}r_{2}})- \theta_{0})^{2}$ This leaves us with $$\frac{\partial \Phi}{\partial r_{1}} = 2k(\arccos(\frac{\mathbf{r_{1}\cdot r_{2}}}{r_{1}r_{2}})- \theta_{0})(-\frac{1}{\sqrt{1-(\frac{\mathbf{r_{1}\cdot r_{2}}}{r_{1}r_{2}})^{2}}})(-\frac{\mathbf{r_{1}\cdot r_{2}}}{r_{1}^{2}r_{2}})$$ $$\frac{\partial \Phi}{\partial r_{2}} = 2k(\arccos(\frac{\mathbf{r_{1}\cdot r_{2}}}{r_{1}r_{2}})- \theta_{0})(-\frac{1}{\sqrt{1-(\frac{\mathbf{r_{1}\cdot r_{2}}}{r_{1}r_{2}})^{2}}})(-\frac{\mathbf{r_{1}\cdot r_{2}}}{r_{1}r_{2}^{2}})$$ The forces that LAMMPs assigns to each atom are then $$\textbf{F}_{1} = (-\frac{\partial \Phi}{\partial r_{1}})(\frac{delx_{1}\hat{x}+dely_{1}\hat{y}+delz_{1}\hat{z}}{r_{1}}) + (-\frac{\partial \Phi}{\partial r_{2}})(\frac{delx_{2}\hat{x}+dely_{2}\hat{y}+delz_{2}\hat{z}}{\frac{\textbf{r}_{1}\cdot\textbf{r}_{2}}{r_{2}}})$$ $$\textbf{F}_{3} = (-\frac{\partial \Phi}{\partial r_{1}})(\frac{delx_{1}\hat{x}+dely_{1}\hat{y}+delz_{1}\hat{z}}{\frac{\textbf{r}_{1}\cdot\textbf{r}_{2}}{r_{1}}}) + (-\frac{\partial \Phi}{\partial r_{2}})(\frac{delx_{2}\hat{x}+dely_{2}\hat{y}+delz_{2}\hat{z}}{r_{2}})$$ $$\textbf{F}_{2} = -\textbf{F}_{1} - \textbf{F}_{3}$$

Any resources, proofs, suggestions on chapters from a textbook, etc. that would help to provide intuition on this problem would be very much appreciated.

EDIT:

I am going to post the derivation of the forces from this bond angle potential as presented by Swope and Ferguson 'Alternative Expressions for Energies and Forces Due to Angle Bending and Torsional Energy' 1992. Hopefully with this information plus the answer provided below by Shern Ren Tee, others will be able to use this post as a starting point for their own calculations.

Consider the following diagram as presented above.

$$\Phi(\theta) = k(\theta - \theta_{0})^{2}$$

We seek the following quantity to derive the forces on atom 1:

$$-\nabla_{1} \Phi(\theta)$$

where $\nabla_{1}$ is the derivative of the $\Phi$ potential holding the position of each other atom (atoms 2 and 3) constant. From here, we use a clever formulation of the chain rule to simplify much of the math for us:

$$-\nabla_{1} \Phi(\theta) = (-\frac{d\Phi}{d\theta})(\frac{d\theta}{d\cos(\theta)})(\nabla_{1}\cos(\theta))$$

The first two terms in this expression are trivial. $$(-\frac{d\Phi}{d\theta}) = 2k(\theta-\theta_{0}) = 2k(\arccos(\frac{\mathbf{r_{1}\cdot r_{2}}}{r_{1}r_{2}})- \theta_{0})$$ $$(\frac{d\theta}{d\cos(\theta)}) = \frac{1}{\frac{\cos(\theta)}{d\theta}} = -\frac{1}{\sin(\theta)} = -\frac{1}{\sqrt{1-\cos^{2}(\theta)}} = (-\frac{1}{\sqrt{1-(\frac{\mathbf{r_{1}\cdot r_{2}}}{r_{1}r_{2}})^{2}}})$$ The last term is slightly trickier. First we must realize the following: $$\nabla_{1} = \hat{x} \frac{\partial}{\partial x_{1}} + \hat{y} \frac{\partial}{\partial y_{1}} + \hat{z} \frac{\partial}{\partial z)}$$ Let us only consider the x component here for the following derivation. We can demonstrate that $$\hat{x}\nabla_{1}\cos(\theta) = \frac{\partial\cos(\theta)}{\partial x_{1}} = \frac{\partial\cos(\theta)}{\partial delx_{1}}\frac{\partial delx_{1}}{\partial x_{1}} + \frac{\partial\cos(\theta)}{\partial dely_{1}}\frac{\partial dely_{1}}{\partial x_{1}} + \frac{\partial\cos(\theta)}{\partial delz_{1}}\frac{\partial delz_{1}}{\partial x_{1}} + \frac{\partial\cos(\theta)}{\partial delx_{2}}\frac{\partial delx_{2}}{\partial x_{1}} + \frac{\partial\cos(\theta)}{\partial dely_{2}}\frac{\partial dely_{2}}{\partial x_{1}} + \frac{\partial\cos(\theta)}{\partial delz_{2}}\frac{\partial delz_{2}}{\partial x_{1}} = \frac{\cos(\theta)}{\partial delx_{1}}\hat{x}$$ Finally, we take that derivative to get the portion of the force in the x direction on atom 1 to see:

$$\frac{\partial\cos(\theta)}{\partial x_{1}} = \frac{\partial}{\partial delx_{1}}(\frac{\mathbf{r}_{1}\cdot\mathbf{r}_{2}}{r_{1}r_{2}}) = -\frac{1}{r_{1}r_{2}}(\frac{\mathbf{r}_{1}\cdot\mathbf{r}_{2}}{r_{1}^{2}}delx_{1}\hat{x} - delx_{2}\hat{x})$$ Or overall we get that the entire $\nabla_{1} \cos(\theta)$ term reduces to $$-\frac{1}{r_{1}r_{2}}(\frac{\mathbf{r}_{1}\cdot\mathbf{r}_{2}}{r_{1}^{2}}\mathbf{r}_{1} - \mathbf{r}_{2})$$ Which is identical to the result obtained from LAMMPs. I skipped over some of the more tedious steps, but I hope that this derivation may help others as a starting point.

Answer 1

Your first question -- writing out analytical forms for forces from potential functions -- is the hardest of the lot. For manual derivations the chain rule is often useful, and your workings already show an application of it. I'm afraid there are no shortcuts there unless you are familiar with a computer mathematics package (like Mathematica or Sage) that can churn out the answers for you.

But proving that (or when) these forces conserve energy and momentum is very easy! Consider the equations of motion (with vector variables as appropriate):

$$ \dot{p} = -\nabla_x U; \hspace{2em} \dot{x} = \frac{p}{m} \tag{1} $$

and the Hamiltonian, or total energy,

$$ H \equiv U + \frac{p^2}{2m}. \tag{2} $$

A quick application of the chain rule gives

$$ \dot{H} = \dot{p} \cdot \nabla_p H + \dot{x} \cdot \nabla_x H \tag{3} $$

and substituting in equation (1) quickly shows that $\dot{H}$ for these equations of motion is zero.

From equation (1), summed over all particles in a system, you can also quickly see that wherever Newton's Third Law holds -- that is, where all forces can be decomposed into action-reaction pairs between particles -- the total force on the system will be zero and momentum will be conserved.

Thus, energy will be conserved in all molecular systems (with thermostatting, more complicated definitions are needed), while momentum often isn't. As a very simple example, two mobile particles connected by a harmonic bond will conserve both energy and momentum; a mobile particle tethered harmonically to a fixed point will conserve energy but not momentum.

Both results are examples of Noether's Theorem, a deep and general result stating that system symmetries result in conserved quantities. For our particular cases, a time-invariant Hamiltonian results in conserved energy and a translation-invariant Hamiltonian results in conserved linear momenta.