Second Derivative of a Rotation Matrix w.r.t. a Rotation vector I'm trying to find the normal mode frequencies of molecular clusters and to do that I require the second derivatives of a Rotation Matrix w.r.t. the components of a rotation vector.
I know how to do this directly using Rodrigue's Rotation Formula, however in my implementation I am parameterizing the rotation matrix with quaternions and I would like to find the second derivatives via quaternions.  
The first derivative is easy to find:
$$ \dfrac{\partial{\textbf{R}}}{\partial{p_k}} = \sum_{i = 0}^{3}\dfrac{\partial{\textbf{R}}}{\partial{q_i}}\dfrac{\partial{q_i}}{\partial{p_k}} $$
where:
$\textbf{R}$ is the 3x3 rotation matrix
$\textbf{q} = (q_0, q_1,q_2, q_3)$ is a unit quaternion
$\textbf{p} = (p_x, p_y,p_z)$ is a rotation vector
$k \in \{x,y,z\} $
$\frac{\partial{\textbf{R}}}{\partial{q_i}}$ is 3x3 matrix and $\frac{\partial{q_i}}{\partial{p_k}}$ is a scalar.
I would like to find $ \dfrac{\partial^{2}{\textbf{R}}}{\partial{p_k}\partial{p_l}}$ 
I've tried using Faa Di Bruno's formula, and this is what I get:
$$ \dfrac{\partial^{2}{\textbf{R}}}{\partial{p_k}\partial{p_l}} = \sum_{i = 0}^{3}\dfrac{\partial^{2}{\textbf{R}}}{\partial{q_i}^{2}}\dfrac{\partial{q_i}}{\partial{p_k}}\dfrac{\partial{q_i}}{\partial{p_l}} + \dfrac{\partial{\textbf{R}}}{\partial{q_i}}\bigg(\dfrac{\partial{}}{\partial{p_k}}\dfrac{\partial{q_i}}{\partial{p_l}}\bigg) $$
I know this is wrong but I just can't figure out how to get to the correct answer.
 A: The answer is:  
Given,
$\textbf{R} \text{ is a 3x3 rotation matrix and can be given as a function of a unit quaternion }\textbf{q} $
$\textbf{q} = (q_0, q_1,q_2, q_3) \text{ is a unit quaternion and can be given as a function of a rotation vector } \textbf{p}$
$\textbf{p} = (p_x, p_y,p_z) \text{ is a rotation vector} $ and $k, l \in \{x,y,z\} $
\begin{align*}
 \dfrac{\partial^{2}{\textbf{R}}}{\partial{p_k}\partial{p_l}} 
&= \dfrac{\partial}{\partial{p_k}}\bigg(\dfrac{\partial{\textbf{R}}}{\partial{p_l}}\bigg) \\
&= \dfrac{\partial}{\partial{p_k}}\bigg(\sum_{i = 0}^{3}\dfrac{\partial{\textbf{R}}}{\partial{q_i}}\dfrac{\partial{{q_i}}}{\partial{p_l}}\bigg) \\
&= \sum_{i = 0}^{3}\dfrac{\partial^{2}{\textbf{R}}}{\partial{p_k}\partial{q_i}}\bigg(\dfrac{\partial{q_i}}{\partial{p_l}}\bigg) + \sum_{i = 0}^{3}\dfrac{\partial{\textbf{R}}}{\partial{q_i}}\bigg(\dfrac{\partial{}}{\partial{p_k}}\dfrac{\partial{q_i}}{\partial{p_l}}\bigg)  \\
&= \sum_{i = 0}^{3}\bigg(\dfrac{\partial{\textbf{R}}}{\partial{p_k}}\bigg)\bigg(\dfrac{\partial}{\partial{q_i}}\bigg)\bigg(\dfrac{\partial{q_i}}{\partial{p_l}}\bigg) + \sum_{i = 0}^{3}\dfrac{\partial{\textbf{R}}}{\partial{q_i}}\bigg(\dfrac{\partial{}}{\partial{p_k}}\dfrac{\partial{q_i}}{\partial{p_l}}\bigg)  \\
&= \sum_{i = 0}^{3}\bigg(\sum_{j = 0}^{3}\dfrac{\partial{\textbf{R}}}{\partial{q_j}}\dfrac{\partial{{q_j}}}{\partial{p_l}}\bigg)\bigg(\dfrac{\partial}{\partial{q_i}}\bigg)\bigg(\dfrac{\partial{q_i}}{\partial{p_l}}\bigg) + \sum_{i = 0}^{3}\dfrac{\partial{\textbf{R}}}{\partial{q_i}}\bigg(\dfrac{\partial{}}{\partial{p_k}}\dfrac{\partial{q_i}}{\partial{p_l}}\bigg)  \\
&= \sum_{i = 0}^{3}\sum_{j = 0}^{3}\dfrac{\partial^{2}{\textbf{R}}}{\partial{q_j}\partial{q_i}}\dfrac{\partial{q_j}}{\partial{p_k}}\dfrac{\partial{q_i}}{\partial{p_l}} + \sum_{i = 0}^{3}\dfrac{\partial{\textbf{R}}}{\partial{q_i}}\bigg(\dfrac{\partial{}}{\partial{p_k}}\dfrac{\partial{q_i}}{\partial{p_l}}\bigg)  
\end{align*}
When I originally posted the question I forgot to take into account that you would need to sum over 2 different quaternions.
