$\rho (\vec{r})$ is the charge density and $\Phi(\vec{r})$ is the electrostatic potential. $$ W= \frac{1}{2} \int \limits_{V} \rho (\vec{r}) \Phi(\vec{r}) dV\quad (i) \\ \vec{\nabla} \cdot \vec{E}(\vec{r})=\frac{\rho(\vec{r})}{\varepsilon_{0}}\quad (ii)$$
Now insert $(ii)$ into $(i)$: $$\frac{\epsilon_{0}}{2} \int (\vec{\nabla} \cdot \vec{E})\Phi d V $$
But in my script the insertion results into the following formula
$$\frac{\epsilon_{0}}{2} \int \Phi \vec{\nabla} \cdot \vec{E} d V =-\frac{\varepsilon_{0}}{2} \int \vec{E} \cdot \vec{\nabla} \Phi d V+\frac{\epsilon_{0}}{2} \int \vec{\nabla} \cdot\left(\vec{E} \Phi\right) d V$$
I know that in $(i)$ $\rho (\vec{r})$ and $\Phi(\vec{r})$ are scalar functions and commutative but is $\Phi (\vec{\nabla} \cdot \vec{E})= (\vec{\nabla} \cdot \vec{E})\Phi$ really the same?
I don't understand the equality: $$\frac{\epsilon_{0}}{2} \int \Phi \vec{\nabla} \cdot \vec{E} d V =-\frac{\varepsilon_{0}}{2} \int \vec{E} \cdot \vec{\nabla} \Phi d V+\frac{\epsilon_{0}}{2} \int \vec{\nabla} \cdot\left(\vec{E} \Phi\right) d V$$
Why can't I just multiply the potential$\Phi$ with the divergence of $\vec{E} $?
By the product rule we have $$ \vec\nabla\cdot(\Phi\vec E) = (\vec\nabla\Phi)\cdot\vec E + \Phi(\vec\nabla\cdot\vec E) . $$
Using this we get $$ W = \frac12 \int \rho(\vec r)\,\Phi(\vec r)\,dV = \frac12 \int \epsilon_0\vec\nabla\cdot\vec E(\vec r)\,\Phi(\vec r)\,dV = \frac{\epsilon_0}{2} \int \vec\nabla\cdot\vec E(\vec r)\,\Phi(\vec r)\,dV \\ = \frac{\epsilon_0}{2} \int \vec\nabla\cdot\left(\vec E(\vec r)\,\Phi(\vec r)\right)\,dV - \frac{\epsilon_0}{2} \int \vec E(\vec r)\cdot\vec\nabla\Phi(\vec r)\,dV . $$