Gradient of an angle in terms of the vertices

Question

Let $\theta(\vec p, \vec q, \vec r)$ be the angle theta between 3D real vectors $(\vec{q}-\vec{p})$ and $(\vec{r} - \vec{p})$. What is a simple expression of $\nabla \theta$ in terms of $\vec{p}$, $\vec{q}$, and $\vec{r}$ that is efficient for computation?

The context is that I'm writing a minimization function that operates on a deformable triangle mesh; I would like to use a harmonic angle potential to prevent significant deformation of the angles in the mesh. The potential energy function for this would be $E(\theta) = K (\theta - \theta_0)^2$ and the gradient in terms of $\theta$ is $\nabla E = 2 K (\theta - \theta_0)$ for a constant $K$ and an initial angle $\theta_0$. Finding the form of $\nabla E$ in terms of the vector values is quite trivial with the exception of finding the form of $\nabla \theta$ in terms of the vectors.

It's easy enough to pop this expression into Mathematica and get a monstrous form of the gradient out in terms of $(x_p, y_p, z_p, \; x_q, y_q, z_q, \; x_r, y_r, z_r)$; if I were only calculating this gradient over a few values this would be fine, but I'll be calculating this for hundreds of thousands of angles, and need the calculation to be very efficient. My intuition is that the derivative in terms of $\vec{q}$ or $\vec{r}$ should be orthogonal to $\vec{q} - \vec{p}$ or $\vec{r} - \vec{p}$, respectively, and that the derivative in terms of $\vec{p}$ should be parallel to the angle bisector. This doesn't particularly help me come up with something more efficient, however. I've found numerous resources that talk about this particular problem (here, for example, section on Angle Potential), usually interms of molecular dynamics simulation, but have been unable to find a simple form of the gradient.

Thanks in advance.

Answer 1

I am also looking answer for the same question. So far I am done with below. I am referring to the equation in Gromacs manual for harmonic angle bond. $$ V(\theta_{ijk})=\frac{1}{2}k_{ijk}^\theta (\theta_{ijk}-\theta_{0})^{2} $$ $$ \theta_{ijk}=cos^{-1} \frac{\overrightarrow{r_{ij}}.\overrightarrow{r_{kj}}}{r_{ij}r_{kj}} $$ $$ F_{i}=-\frac{\partial V(\theta)}{\partial r_{i}} $$ $$ F_{i}= -k_{ijk}^\theta.(\theta_{ijk}-\theta_{0}).\frac{\partial \theta_{ijk}}{\partial r_i} $$ $$ \frac{\partial \theta_{ijk}}{\partial r_i}=\frac{-1}{\sqrt{1-u^{2}}}.\frac{\partial u}{\partial r_i} $$ where $$u=\frac{\overrightarrow{r_{ij}}.\overrightarrow{r_{kj}}}{r_{ij}r_{kj}} $$ Now, $$ \frac{\partial u}{\partial r_i}=\frac{1}{r_{ij}r_{kj}}(\overrightarrow{r_{ij}}.\frac{d\overrightarrow{r_{kj}}}{dr_i}+\frac{d\overrightarrow{r_{ij}}}{dr_i}.\overrightarrow{r_{kj}} ) ?? $$ I am stuck in the last part...Anybody with better grip in maths can please point out the missing detail...

Answer 2

This website gives an answer to the problem; see equations 4.108 - 4.110: http://www.mbnexplorer.com/documentation/3-energy-and-force-calculation/33-molecular-mechanics-potential

In case the link disappears, here are the equations, translated into the notation of the original question; note that these give the gradient, as the question asked, while the above website, which is about the molecular dynamics problem in particular, gives the force ($-\nabla$).

$$ E(\vec{p}, \vec{q}, \vec{r}) = K\left( \theta(\vec{p}, \vec{q}, \vec{r}) - \theta_0\right)^2\\ \nabla E(\vec{p}, \vec{q}, \vec{r}) = \left( \frac{\partial E(\vec{p}, \vec{q}, \vec{r})}{\partial \vec{p}}, \frac{\partial E(\vec{p}, \vec{q}, \vec{r})}{\partial \vec{q}}, \frac{\partial E(\vec{p}, \vec{q}, \vec{r})}{\partial \vec{r}}\right)\\ \frac{\partial E(\vec{p}, \vec{q}, \vec{r})}{\partial \vec{p}} = -\frac{\partial E(\vec{p}, \vec{q}, \vec{r})}{\partial \vec{q}} - \frac{\partial E(\vec{p}, \vec{q}, \vec{r})}{\partial \vec{r}}\\ \frac{\partial E(\vec{p}, \vec{q}, \vec{r})}{\partial \vec{q}} = \frac{2K\left( \theta(\vec{p},\vec{q},\vec{r}) - \theta_0\right)}{||\vec{q}-\vec{p}||\,\sin\left(\theta(\vec{p},\vec{q},\vec{r})\right)} \left( \frac{\vec{r}-\vec{p}}{||\vec{r}-\vec{p}||} - \frac{\vec{q}-\vec{p}}{||\vec{q}-\vec{p}||}\cos\left(\theta(\vec{p},\vec{q},\vec{r})\right) \right)\\ \frac{\partial E(\vec{p}, \vec{q}, \vec{r})}{\partial \vec{r}} = \frac{2K\left( \theta(\vec{p},\vec{q},\vec{r}) - \theta_0\right)}{||\vec{r}-\vec{p}||\,\sin\left(\theta(\vec{p},\vec{q},\vec{r})\right)} \left( \frac{\vec{q}-\vec{p}}{||\vec{q}-\vec{p}||} - \frac{\vec{r}-\vec{p}}{||\vec{r}-\vec{p}||}\cos\left(\theta(\vec{p},\vec{q},\vec{r})\right) \right) $$

While I haven't derived the correct solution, I did plop the solution into a Mathematica instance and verified that it agreed with the numerical solution for a wide range of values.

Answer 3

This is a very late answer, but I just leave it here for others, since it seems these things are not readily available online:

$$ \theta = acos (\frac{\vec r_{ij} \cdot \vec r_{kj}}{r_{ij} r_{kj}}) $$

Let's look at the derivative with respect to $\vec r_i$:

$$ \frac{\partial \theta}{\partial \vec r_i} = \frac{\partial \theta}{\partial r_{ij}} \frac{\partial r_{ij}}{\partial \vec r_{ij}} \frac{\partial \vec r_{ij}}{\partial \vec r_{i}} = \frac{-1}{\sqrt{1-(\frac{\vec r_{ij} \cdot \vec r_{kj}}{r_{ij} r_{kj}})^2}} (\frac{\partial \vec r_{ij}}{\partial r_{ij}} \cdot \frac{\vec r_{kj}}{r_{ij} r_{kj}} + \frac{-1}{r_{ij}^2} \frac{\vec r_{ij} \cdot \vec r_{kj}}{r_{kj}}) \frac{\partial r_{ij}}{\partial \vec r_{ij}} \frac{\partial \vec r_{ij}}{\partial \vec r_{i}} $$

Since $\vec r_{ij} = \vec r_i - \vec r_j$, the last derivative term is simply identity. Also, for $\vec r = (x,y,z)$, $r = \sqrt{x^2+y^2+z^2}$, so the second last derivative is:

$$ \frac{\partial r}{\partial \vec r} = \frac{1}{2 \sqrt{x^2+y^2+z^2}}(2x, 2y, 2z) = \frac{\vec r}{r} $$

Thus: $$ \frac{\partial \theta}{\partial \vec r_i} = \frac{- r_{ij} r_{kj}}{\sqrt{r_{ij}^2 r_{kj}^2 - (\vec r_{ij} \cdot \vec r_{kj})^2}} (\frac{\vec r_{kj}}{r_{ij} r_{kj}} - \frac{\vec r_{ij} \cdot \vec r_{kj}}{r_{ij}^2 r_{kj}} \frac{\vec r_{ij}}{r_{ij}}) = \frac {-1}{r_{ij}} \frac{r_{ij}^2 \vec r_{kj} - (\vec r_{ij} \cdot \vec r_{kj}) \vec r_{ij}} {\sqrt{r_{ij}^4 r_{kj}^2 - (\vec r_{ij} \cdot \vec r_{kj})^2 r_{ij}^2}} $$

The numerator is the vector triple product, so the final result is: $$ \frac{norm(\vec r_{ij} \times (\vec r_{ij} \times \vec r_{kj}))}{r_{ij}} $$

You can apply the same idea to find the derivative with respect to $\vec r_k$ and $\vec r_j$. Also, this also works with dihedral angle formula, just with more things to take care of.