As you suspect, the mechanical energy equation you derive is missing terms that account for viscous effects. The other terms you derive are correct.
The RHS should be $\nabla \cdot \mathbf{\sigma}$, the divergence of the stress tensor, and you have edited accordingly.
This clearly is the correct term to include because the rate of increase of momentum in a fixed volume $V$ of fluid should be related to the viscous stress force of the surrounding fluid acting on the bounding surface of the volume $\partial V$,
$$\text{rate of momentum change } = \int_{\partial V}\mathbf{\sigma} \cdot \mathbf{n} \, dS = \int_V \nabla \cdot \mathbf{\sigma} \, dV.$$
Your error stems from incorrectly forming the dot product of the velocity and the viscous term $\eta\,\nabla \cdot [\nabla \mathbf{v} + \nabla \mathbf{v}^T] = \eta \, \nabla^2 \mathbf{v} + \eta \, \nabla (\nabla \cdot \mathbf{v})$.
In particular, you are missing $ \mathbf{v} \cdot [\eta\,\nabla \cdot \nabla \mathbf{v}]$ which should give rise to the viscous dissipation term, representing the action of viscous forces in converting useful mechanical energy into heat.
The correct treatment of that term is
$$\tag{*} \eta\,\mathbf{v} \cdot\,\nabla \cdot [\nabla \mathbf{v} + \nabla \mathbf{v}^T] = \underbrace{-\eta \nabla \mathbf{v} \,\mathbf{:}\, \nabla \mathbf{v}}_{\text{viscous dissipation}} + \frac{\eta}{2}\nabla^2(|\mathbf{v}|^2) + \eta \, \mathbf{v} \cdot \nabla (\nabla \cdot \mathbf{v}), $$
where the double dot product represents $\sum_{i=1}^3 \sum_{j=1}^3 |\frac{\partial v_i}{\partial x_j}|^2$, a norm of the velocity gradient, and $|\mathbf{v}|$ is the usual Euclidean norm of the velocity vector.
For an easy derivation, we can expand $\mathbf{v} \cdot\,\nabla \cdot [\nabla \mathbf{v}+ \nabla \mathbf{v}^T] $ using Cartesian coordinates and the Einstein convention (repeated index implies summation over that index). First we obtain
$$\nabla \cdot [\nabla \mathbf{v}+ \nabla \mathbf{v}^T]_{i} = \partial_j(\partial_j v_i + \partial_iv _j) = \partial_j(\partial_j v_i) + \partial_i(\partial_jv_j),$$
and forming the dot product we get
$$ \mathbf{v} \cdot\,\nabla \cdot [\nabla \mathbf{v} + \nabla \mathbf{v}^T] = v_i\,\partial_j(\partial_j v_i) + v_i \,\partial_i(\partial_jv_j) \\ = -\partial_jv_i \, \partial_j v_i + \partial_j[v_i(\partial_j v_i)] + v_i \,\partial_i(\partial_jv_j) \\ = -\partial_jv_i \, \partial_j v_i + \frac{1}{2}\partial_j\,\partial_j (v_iv_i) + v_i \,\partial_i(\partial_jv_j), $$
which is the component form of the coordinate-free term on the RHS of (*).
Best Answer
The definition below may not be appropriate, but defining $\mathbf{\vec{r}}$ as:
$$\mathbf{\vec{r}}_{ij} = x \mathbf{\hat{i}} + y \mathbf{\hat{j}} + z \mathbf{\hat{k}}$$
Where $\mathbf{\vec{r}}_{ij}$ defines the vector between atoms $i, j$, and the values $\mathbf{\hat{i}}, \mathbf{\hat{j}}, \mathbf{\hat{k}}$ are unit vectors.
Gives you:
$$r_{ij} = |\mathbf{\vec{r}}_{ij}| = \sqrt{x^2+y^2+z^2} = \left(x^2+y^2+z^2\right)^{1/2}$$
Then defining $\nabla$ as:
$$\nabla = \frac{\partial}{\partial x} \mathbf{\hat{i}} + \frac{\partial}{\partial y} \mathbf{\hat{j}} + \frac{\partial}{\partial z} \mathbf{\hat{k}} $$
Gives you:
$$\nabla u(r_{ij}) = \frac{\partial u(r_{ij})}{\partial x} \mathbf{\hat{i}} + \frac{\partial u(r_{ij})}{\partial y} \mathbf{\hat{j}} + \frac{\partial u(r_{ij})}{\partial z} \mathbf{\hat{k}} $$
This allows you to use the chain rule for each partial derivative, so for each partial derivative we have a variant of what's shown below for the partial with respect to $x$:
$$\frac{\partial u(r_{ij})}{\partial x} = \frac{\partial r_{ij}}{\partial x}\frac{\partial u(r_{ij})}{\partial r_{ij}}$$
Letting $r_{ij} = r$ for notational simplicity we have:
$$\frac{\partial u(r)}{\partial x} = \frac{\partial r}{\partial x}\frac{\partial u(r)}{\partial r}$$
The first part of the right hand side of the equation can be calculated for $x$ as:
$$\begin{eqnarray}\frac{\partial r}{\partial x} = \frac{\partial \left(x^2 + y^2 + z^2\right)^{1/2}}{\partial x} &=& \left(2x\right) \left(\frac{1}{2}\right) \left(x^2 + y^2 + z^2\right)^{-1/2} \\ \\ &=& \frac{x}{\left(x^2 + y^2 + z^2 \right)^{1/2}}\end{eqnarray}$$
As $r = \left(x^2+y^2+z^2\right)^{1/2}$ this simplifies to:
$$\frac{\partial r}{\partial x} = \frac{x}{r}$$
Similar equations for the partial with respect to $y$ and $z$ can be found. This leaves the second part of the right hand side equation:
$$u(r) = D_0\left[ e^{-2\alpha(r-r_0)} -2e^{-\alpha(r-r_0)}\right]$$
Where again we have let $r_{ij} = r$, taking the partial with respect to r we have:
$$\frac{\partial u(r)}{\partial r} = D_0\left[ \left(-2\alpha\right)e^{-2\alpha(r-r_0)} - (-2\alpha)e^{-\alpha(r-r_0)}\right]$$
As $\frac{\partial u(r)}{\partial r}$ will be the same for all partials with respect to $x, y, z$ it can be factored out of the final equation and we have:
$$\nabla u(r) = \frac{\partial u(r)}{\partial r}\left(\frac{\partial r}{\partial x}\mathbf{\hat{i}} + \frac{\partial r}{\partial y}\mathbf{\hat{j}} + \frac{\partial r}{\partial z}\mathbf{\hat{k}}\right)$$
Subbing in the partials with respect to $x, y, z$ along with $\frac{\partial u(r)}{\partial r}$ you get:
$$\begin{eqnarray}\nabla u(r) &=& D_0\left[ \left(-2\alpha\right)e^{-2\alpha(r-r_0)} - (-2\alpha)e^{-\alpha(r-r_0)}\right]\left(\frac{x}{r}\mathbf{\hat{i}} + \frac{y}{r}\mathbf{\hat{j}} + \frac{z}{r}\mathbf{\hat{k}}\right) \\ \\ &=& \frac{-\left(2 \alpha\right)D_0}{r} \left[ e^{-2\alpha(r-r_0)} - e^{-\alpha(r-r_0)}\right]\left(x\mathbf{\hat{i}} + y\mathbf{\hat{j}} + z\mathbf{\hat{k}}\right)\end{eqnarray}$$
As $\mathbf{\vec{r}}_{ij} = x \mathbf{\hat{i}} + y \mathbf{\hat{j}} + z \mathbf{\hat{k}}$ you have:
$$\nabla u(r) = \frac{-\left(2 \alpha\right)D_0}{r} \left[ e^{-2\alpha(r-r_0)} - e^{-\alpha(r-r_0)}\right]\mathbf{\vec{r}}$$
Which finally simplifies to:
$$\nabla u(r) = -\left(2 \alpha\right)D_0 \left[ e^{-2\alpha(r-r_0)} - e^{-\alpha(r-r_0)}\right]\mathbf{\hat{r}}$$
Where $\mathbf{\hat{r}} = \frac{\mathbf{\vec{r}}}{r}$ and is a unit vector in the direction of $\mathbf{\vec{r}}$, here the notation has also been simplified and $\mathbf{\vec{r}} = \mathbf{\vec{r}}_{ij}$ from the original definition.