1) The first derivation is correct. When you make the second one for $\dot x$ written as $\dot x = \dfrac{d x}{d t}$ you make an error when using full derivatives instead of partial ones. In the latter case you would have obtained:
$$
\dot x=\dfrac{d x}{d t}=\dfrac{\partial x}{\partial \theta}\dfrac{d \theta}{d t} +\dfrac{\partial x}{\partial \phi}\dfrac{d \phi}{d t}+\dfrac{\partial x}{\partial r}\dfrac{d r}{d t}
$$
Then, as $x=r \cos \theta$, you can insert $\dfrac{\partial x}{\partial \phi}=0, \dfrac{\partial x}{\partial \theta}=-r \sin \theta$, $\dfrac{\partial x}{\partial r}=\cos \theta$ into the last formula to obtain the same result you had before.
2) The derivatives you are using here are not covariant derivatives, they are always partial. The derivative $\dfrac{d}{d t}$ is expressed through partial derivatives, in contrast to full derivative $\dfrac{D}{D t}$, which is expressed through covariant derivatives. As you are actually not using covariant derivatives in your derivation, you don't have to worry about the basis vectors.
3) You could make you derivation alternatively by rewriting your original equations in a covariant form $M\dfrac{D^2 x^i}{D t^2}+C\dfrac{D x^i}{D t}+Kx^i=0$, where $x^i$ stands for $(x,y,z)$, stating that the equations of motion do not depend on the choice of coordinates, rewriting the equation in spherical coordinates as $M\dfrac{D^2 \tilde x^i}{D t^2}+C\dfrac{D \tilde x^i}{D t}+K \tilde x^i=0$, where $\tilde x = (r, \theta, \phi)$, and expanding all the derivatives $\dfrac{D}{D t}$ through partial derivatives and Christoffel symbols in spherical coordinates. You would have to take care of the basis, as the result would be in natural basis, whereas the equations of motion are usually written in orthonormal basis.
4) Another powerful way to transform coordinates would be by using lagrangian form of the equations of motion. This would be possible for the case of conservative systems only (no damping terms). In this case you could just rewrite the lagrangian of your system in spherical coordinates and write lagrangian equations from it in a more or less straightforward manner.
Best Answer
The formula $$ \sum_{i=1}^3 p_i q_i $$ for the dot product obviously holds for the Cartesian form of the vectors only. The proposed sum of the three products of components isn't even dimensionally correct – the radial coordinates are dimensionful while the angles are dimensionless, so they just can't be added.
One can't simplify the calculation much. At most, one may realize that the inner product will only depend on $\phi_1,\phi_2$ through their difference $\phi_1-\phi_2$ because one may use the rotational symmetry around the $z$-axis to set e.g. $\phi_2$=0.
While doing so, we may set $\phi_2=0$ i.e. $y_2=0$ and the inner product reduces to $$ x_1 x_2 + z_1 z_2 = r_1r_2 (\sin\theta_1\sin\theta_2 \cos\phi_1 + \cos\theta_1\cos\theta_2) $$ We may restore the form for a general rotation by replacing $\phi_1$ in the formula above by $\phi_1-\phi_2$ to get the inner product $$ r_1r_2 (\sin\theta_1\sin\theta_2 \cos(\phi_1-\phi_2) + \cos\theta_1\cos\theta_2) $$ which is the same as your formula because $\cos(a-b)=\cos a\cos b +\sin a \sin b$.