Obtaining Brans-Dicke theory scalar (wave) equation I have trouble with obtaining d'Alambert equation for scalar field in Brans-Dicke gravity (http://www.scholarpedia.org/article/Jordan-Brans-Dicke_Theory).
B-D gravity langrangian density is given by:
$$\mathcal{L}=\phi R + 16\pi\mathcal{L}_{matter} -\frac{\omega(\phi_{,\alpha}\phi^{,\alpha}
)}{\phi},$$
with corresponding variational principle:
$$\delta\int (-g)^{1/2}[\phi R + 16\pi\mathcal{L}_{matter} -\frac{\omega(\phi_{,\alpha}\phi^{,\alpha}
)}{\phi}]d^4x = 0.$$
When i'm using standard Euler-Langrange equation for $\phi$:
$\frac {\partial\mathcal{L}}{\partial\phi}-\partial_\mu(\frac{\partial\mathcal{L}}{\partial(\partial_\mu\phi)})=0$
I get:
$\frac {\partial\mathcal{L}}{\partial\phi}=R + (-\omega\frac {\partial{\frac{\phi_{,\alpha}\phi^{,\alpha}}{\phi}}}{\partial\phi})=R-(-\frac{\omega(\phi_{,\alpha}\phi^{,\alpha}
)}{\phi^2})=R + \frac{\omega(\phi_{,\alpha}\phi^{,\alpha}
)}{\phi^2}$ 
and:
$\partial_\mu(\frac{\partial\mathcal{L}}{\partial(\partial_\mu\phi)})=-\frac{\omega}{\phi} \partial_\mu(\frac{\partial\mathcal{\frac{\phi_{,\alpha}\phi^{,\alpha}}{\phi}}}{\partial(\partial_\mu\phi)})=-\frac{\omega}{\phi}  g^{\alpha\beta}(\partial_\mu(\frac{\partial_\alpha\phi\partial_\beta\phi}{\partial(\partial_\mu\phi)}))=-\frac{\omega}{\phi} g^{\alpha\beta}\partial_\mu(\delta^\mu_\alpha\partial_\beta\phi+\delta^\mu_\beta\partial_\alpha\phi)=-\omega \partial_\mu[g^{\mu\beta}\partial_\beta\phi+g^{\alpha\mu}\partial_\alpha\phi]=-\frac{\omega}{\phi} \partial_\mu[g^{\beta\mu}\partial_\beta\phi+g^{\alpha\mu}\partial_\alpha\phi]
=-\frac{\omega}{\phi}  \partial_\mu2g^{\mu\beta}\partial_\beta\phi=-2\frac{\omega}{\phi} g^{-1/2} \partial_\mu[g^{1/2}g^{\mu\beta}\partial_\beta\phi]=-2\frac{\omega}{\phi} \Box\phi$
Which after plugging into E-L equation gives:
$0=R+\omega\frac{\phi_{,\alpha}\phi^{,\alpha}}{\phi^2}+(2\omega\Box\phi)/\phi$
After multiplying by $\phi/2\omega$ i arrived at:
$0=\frac{\phi R}{2\omega}+\frac{\phi_{,\alpha}\phi^{,\alpha}}{2\phi}+\Box\phi$
However correct equation should be (from article linked above):
$0=\frac{\phi R}{2\omega}-\frac{\phi_{,\alpha}\phi^{,\alpha}}{2\phi}+\Box\phi$
How to do it correctly?
 A: There is a mistake in the beginning, you should use E-L equations in curved spacetime. Here is a detailed calculation
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] $$
EoM for BD scalar field
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}
 Let us evaluate the $\phi$ differentiation,
 \begin{align}
 \frac{\partial\mathcal{L}}{\partial\phi} = R + \frac{\omega}{\phi^2} g^{\mu \nu} \nabla_{\mu} \phi \nabla_\nu \phi - \frac{\partial V}{\partial \phi}
 \end{align}
 and evaluating RHS yields, 
 \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}
