Writing down the derivative of the total energy $E$ with respect to the position $\vec{R}_i$ of ion $i$
\begin{equation}
\frac{dE}{d\vec{R_i}}=
- \sum_j \Phi_j^* \frac{\partial H}{\partial \vec{R}_i} \Phi_j
-+\sum_j \frac{\partial \Phi_j^*}{\partial \vec{R}_i} H \Phi_j
-+\sum_j \Phi_j^* H \frac{\partial \Phi_j}{\partial \vec{R}_i}
+ \sum_j \langle \Phi_j | \frac{\partial H}{\partial\vec{R}_i} | \Phi_j \rangle
++\sum_j \langle \frac{\partial \Phi_j}{\partial\vec{R}_i} | H \Phi_j \rangle
++\sum_j \langle \Phi_j H | \frac{\partial \Phi_j}{\partial \vec{R}_i} \rangle
\text{ ,}
\end{equation}
indeed reveals a contribution to the change in total energy due to the change of the wave functions $\Phi_j$.