EoM for scalar field in Brans-Dicke Theory

The action is given by $$S^{(BD)} = \int d^4 x \sqrt{|g|} \left[ \phi R - \frac{\omega}{\phi} g^{\mu \nu} \, \nabla_\mu \phi \nabla_\nu \phi - V(\phi) \right]$$

I am trying to vary with respect to $$\phi$$ using Euler - Lagrange equations in curved spacetime, i.e. $$\underbrace{\frac{\mathcal{\partial L}}{\partial \phi}}_{I} - \underbrace{\nabla_\mu \left( \frac{\mathcal{\partial L}}{\partial (\nabla_\mu \phi)} \right)}_{II} = 0$$

I obtained the following $$I: R + \frac{\omega}{\phi^2}g^{\mu \nu} \, \nabla_\mu \phi \nabla_\nu \phi - \frac{\partial V}{\partial \phi}$$

\begin{align} \begin{aligned} II&: \nabla_\mu \left[ - \frac{\omega}{\phi} g^{\mu \nu} \left( \nabla_\nu \phi + \underbrace{\nabla_\mu \phi \, \delta^\mu{}_\nu}_{\nabla_\nu \phi}\right)\right] \\[3 ex] &= - \nabla_\mu \left( \frac{2\omega}{\phi} \underbrace{g^{\mu \nu} \, \nabla_\nu \phi}_{\nabla^\mu \phi}\right) \\[3 ex] &= -2 \omega \, \nabla_\mu \left( \frac{1}{\phi} \nabla^{\mu} \phi\right) \\[3 ex] &= -2 \omega \left( \nabla_\mu \frac{1}{\phi} \cdot \nabla^\mu \phi + \frac{1}{\phi} \nabla_\mu \nabla^\mu \phi\right)\\[3 ex] & = -2 \omega \left( \nabla_\mu \frac{1}{\phi} \cdot \nabla^\mu \phi + \frac{1}{\phi} \Box \phi\right) \end{aligned} \end{align} According to literature, EoM should have been \begin{align} I - II = R - \frac{\omega}{\phi^2}g^{\mu \nu} \, \nabla_\mu \phi \nabla_\nu \phi - \frac{\partial V}{\partial \phi} + \frac{2\omega}{\phi} \Box \phi \end{align}

So I have an extra term $$\frac{2\omega}{\phi} \nabla_\mu \frac{1}{\phi} \cdot \nabla^\mu \phi$$

I've tried to write step-by-step, any help would be appreciated

Euler - Lagrange Equations in curved spacetime \begin{align} \frac{\partial\mathcal{L}}{\partial\phi}=\frac{1}{\sqrt{-g}}\partial_{\mu}\left[\sqrt{-g}\frac{\partial\mathcal{L}}{\partial\left(\partial_{\mu}\phi\right)}\right] \end{align} OR \begin{align} \frac{\partial\mathcal{L}}{\partial\phi}=\nabla_{\mu}\left[\frac{\partial\mathcal{L}}{\partial\left(\partial_{\mu}\phi\right)}\right] \end{align}
\begin{align} \frac{\partial\mathcal{L}}{\partial\phi} = \phi R + \frac{\omega}{\phi^2} g^{\mu \nu} \nabla_{\mu} \phi \nabla_\nu \phi - \frac{\partial V}{\partial \phi} \end{align}
\begin{align} \begin{aligned} \nabla_{\mu}\left[\frac{\partial\mathcal{L}}{\partial\left(\partial_{\mu}\phi\right)}\right] &= \nabla_{\mu} \left[ - \frac{\omega}{\phi} g^{\mu \nu} \left(\nabla_\nu \phi + \underbrace{\nabla_{\mu} \phi \, \delta^\mu{}_\nu}_{\nabla_\nu \phi} \right)\right] \\[1ex] &= \nabla_{\mu} \left[ - \frac{\omega}{\phi} g^{\mu \nu} \left(2\nabla_\nu \phi \right)\right]\\[1ex] &= -2 \omega \nabla_{\mu} \left( \frac{1}{\phi} g^{\mu \nu} \nabla_{\nu} \phi \right) \\[1ex] &= -2 \omega \nabla_{\mu} \left( \frac{1}{\phi} \nabla^{\mu} \phi \right) \\[1ex] &= -2 \omega \left[ \left(\nabla_{\mu} \frac{1}{\phi}\right) \nabla^\mu \phi + \frac{1}{\phi} \underbrace{\nabla_\mu \nabla^\mu \phi}_{\Box \phi} \right] \\[1ex] &= -2 \omega \left[ \left(-\frac{1}{\phi^2} \nabla_{\mu} \phi \right) \nabla^\mu \phi + \frac{1}{\phi} \Box \phi \right]\\[1ex] &= \frac{2 \omega}{\phi^2} g^{\mu \nu} \nabla_{\mu} \phi \, \nabla_{\nu} \phi - \frac{2\omega}{\phi} \Box \phi \end{aligned} \end{align}
Therefore, \begin{align} \begin{aligned} \frac{\partial\mathcal{L}}{\partial\phi} - \nabla_{\mu}\left[\frac{\partial\mathcal{L}}{\partial\left(\partial_{\mu}\phi\right)}\right] = R + \frac{\omega}{\phi^2} g^{\mu \nu} \nabla_{\mu} \phi \, \nabla_{\nu} \phi - \frac{\partial V}{\partial \phi} + \frac{2 \omega}{\phi} \Box \phi \end{aligned} \end{align} \begin{align} \boxed{R - \frac{\omega}{\phi^2} g^{\mu \nu} \nabla_{\mu} \phi \, \nabla_{\nu} \phi - \frac{\partial V}{\partial \phi} + \frac{2 \omega}{\phi} \Box \phi = 0} \end{align}