Stuck on proving Bochner's formula 
Let $u$ be a smooth function over a Riemannian manifold $(M,g)$. Prove
the following Bochner's formula
\begin{align*}
  \frac{1}{2} \Delta |\nabla u|^2= |\nabla^2 u|+\text{Ric}(\nabla u,\nabla u)+g(\nabla\Delta u,\nabla u)
 \end{align*}
where $\Delta$ is the Laplacian and $|\cdot|^2=g(\cdot,\cdot)$.

My attempt
It suffices to prove the formula at any point $p$.
Choose a normal frame $\{e_i\}_{i=1}^n$ near $p$. Then for any smooth function $f$ we have
\begin{align*}
  \Delta f
  =(\nabla^2f)(e_i,e_i)
  =(\nabla_{e_i} df)(e_i)
  = e_i (df,e_i)- (\nabla_{e_i}e_i,df)
 \end{align*}
and hence
\begin{align*}
  \Delta f(p)=(e_ie_if)(p)
 \end{align*}
Thus the first term is
\begin{align*}
  \boxed{\frac{1}{2} \Delta |\nabla u|^2(p)
  = \frac{1}{2} (e_ie_i (e_ju)^2)(p) 
  = (e_ie_j u)^2(p)+e_ju(p)\cdot e_ie_ie_ju(p)}\tag{1}
 \end{align*}
Now we compute the second term. Note that
\begin{align*}
  |\nabla^2 u|
  = \left(  (\nabla^2 u)(e_i,e_j)  \right)^2
  = \left(  (\nabla_{e_j} du)(e_i)  \right)^2
  = \left(  e_j (du,e_i) - (du,\nabla_{e_j}e_i) \right)^2
 \end{align*}
and hence
\begin{align*}
  \boxed{|\nabla^2 u|(p)= (e_je_i u)^2(p)}\tag{2}
 \end{align*}
Now we compute the third term. Note that
\begin{align*}
  \text{Ric}(\nabla u,\nabla u)
  &= R(e_i,\nabla u,e_i,\nabla u)
  = R(e_i,e_j,e_i,e_k) e_ju \cdot e_k u\\
  &= \left(\left<\nabla_{[e_i,e_j]} e_i,e_k\right>-\left<[\nabla_{e_i},\nabla_{e_j}]e_i,e_k\right>\right)
  e_ju \cdot e_k u
 \end{align*}
and hence
\begin{align*}
  \boxed{\text{Ric}(\nabla u,\nabla u)(p)
  =\left<\nabla_{e_j}\nabla_{e_i}e_i,e_k\right>e_ju \cdot e_k u\Big|_p
   -\left<\nabla_{e_i}\nabla_{e_j}e_i,e_k\right>e_ju \cdot e_k u\Big|_p}\tag{3}
 \end{align*}
Now we compute the last term. Note that
\begin{align*}
  g(\nabla\Delta u,\nabla u)
  &= \left<\nabla\Delta u,e_i\right>\left<\nabla u,e_i\right>\\
  &= e_j \left(e_i e_i u- \left<\nabla_{e_i}e_i,\nabla u\right>\right) e_ju\\
  &= e_j e_i e_i u\cdot e_ju- e_j\left<\nabla_{e_i}e_i,\nabla u\right>\cdot e_ju
 \end{align*}
and hence
\begin{align*}
  \boxed{
  g(\nabla\Delta u,\nabla u)(p)
  =(e_j e_i e_i u)(p)\cdot e_ju(p)
  -\left<\nabla_{e_j}\nabla_{e_i}e_i,e_k \right>\Big|_pe_ku(p)\cdot e_ju(p)}\tag{4}
 \end{align*}
Thus $(1)-(2)-(3)-(4)=0$ is equivalent to that
$$
\left<\nabla_{e_i}\nabla_{e_j}e_i,e_k\right>e_ju \cdot e_k u\Big|_p=0
$$
I got stuck on proving this. What's wrong with my proof? Any helps will be highly appreciated!
 A: The problem is, $e_ie_ie_j u$ and $e_je_ie_i u$ do not cancel. In fact, at $p$ we have
$e_ie_ie_j u-e_je_ie_i u=e_ie_ie_j u-e_ie_je_i u=e_i([e_i, e_j]u)$; now although $[e_i, e_j]=0$ at $p$, its derivative may not vanish; so $ e_i([e_i, e_j]u)$ may not be zero.
Now we compute
$(e_ie_ie_j u-e_je_ie_i u)(e_ju)=\Big(e_i([e_i, e_j]u)\Big)e_ju$; this is
$$e_i\langle [e_i, e_j], \nabla u\rangle e_ju 
=e_i\langle \nabla_{e_i}e_j-\nabla_{e_j}e_i, \nabla u\rangle  \cdot e_ju
=\langle \nabla_{e_i}\nabla_{e_i}e_j-\nabla_{e_i}\nabla_{e_j}e_i, \nabla u\rangle \cdot e_j u;
$$
the second term $-\langle \nabla_{e_i}\nabla_{e_j}e_i, \nabla u\rangle \cdot e_j u$ cancels the last remaining term you wrote; as for the first term, an easy way to get rid of it is, assume the frame is constructed like this: pick an orthonormal frame $e_i$ at $p$, then along any geodesic $\gamma$ from $p$ we parallel translate this frame along $\gamma$; this gives
a frame field locally defined near $p$. In particular, $\nabla_{e_i}e_j=0$ along the geodesic $\gamma_i$ at $p$ with initial tangent vector  $e_i$, therefore $\nabla_{e_i}\nabla_{e_i} e_j=0$ along this geodesic $\gamma_i$, therefore at $p$ we have
$\nabla_{e_i}\nabla_{e_i} e_j=0$.
