Energy in dynamical variational principle In quantum mechanics we use variational principle in order to find approximate expression for the ground state. Lets assume our probe wavefunction $|\Psi\rangle$ can be expanded in orthonormal basis 
$$|\Psi\rangle = \sum\limits_{n}f_n |n\rangle$$
Variational ansatz dictates minimization of the energy functional
$$E[f_n,f_n^*] = \langle\Psi|\hat{H}|\Psi\rangle$$
with constraint $\langle \Psi|\Psi\rangle=1$ ($f_n^*$ is a complex conjugate). Taking derivative we have equations for the coefficients:
$$\frac{\partial E}{\partial f_n^*} = 0.$$
There is also dynamical variational principle where one minimizes Schrodinger action
$$S = \int dt \mathcal{L}$$
where
$$\mathcal{L} = \langle\Psi(t)|i\hbar\partial_t - \hat{H}|\Psi(t)\rangle$$
Using Euler-Lagrange equations we get differential equations for $f_{n}$:
$$\frac{\partial \mathcal{L}}{\partial f_n^*} - \frac{d}{dt}\frac{\partial\mathcal{L}}{\partial \dot{f}_n^*}=0$$
My question is whether or not energy $E(t) = \langle\Psi(t)|\hat{H}|\Psi(t)\rangle$ defined with coefficients $f_n(t)$ derived from Euler-Lagrange equations is a conserved quantity i.e. $d E(t)/dt = 0$?
What I have in mind is the Bose-Hubbard Hamiltonian
$$\hat{H} = -J\sum\limits_{<i,j>}\hat{a}_i^{\dagger}\hat{a}_j + \frac{U}{2}\sum\limits_i\hat{n}_i(\hat{n}_i-1)-\mu\sum\limits_{i}\hat{n}_i$$
with variational ansatz:
$$|\Psi\rangle = \bigotimes\limits_{i}|\psi_i\rangle,\ \ |\psi_i\rangle = \sum\limits_{n=0}^{n_F}f_{n}^{(i)}|n\rangle_{i}$$
In this case Coefficients $f_n$ are not a solution of the Schrodinger equation $i\hbar \partial_t |\Psi\rangle = \hat{H} |\Psi\rangle$.
 A: If you take the equation of motion on the Schrodinger Field Action, you are simple imposing this equation:
$$
\frac{\delta S}{\delta f_{n}}|_{Bulk}=0\,\,\,\,\rightarrow \,\,\,\,\, H_{n}^{m} f_m=i\hbar\frac{df_n (t)}{dt}
$$
Note that I used the word Bulk below the functional derivative. This is because, functional derivative acts only in the interior of the functional, and don't tells you how the boundary terms behaves. In this problem the boundary term are related to the initial state of the system. This is what you want to find. The initial state that you want to find is precisely the state that conserves the probabilities:
$$
\frac{d(f^{*}_{n}f^{n})}{dt}(t)=0
$$
So, the answer that the variation inside the Bulk of the action give us is that if you want such state you need try to diagonalize $H_{mn}$. Very trivial answer. Actually what you want more is the lowest eigenvalue of $H$, so you already know that. The answer for your question is that the equation of motion is not necessary to specify your state. The equation of motion only tells you how some state evolutes through time.
If you are interested in functional applications to vacuum state you can use the idea that in the limit $T\rightarrow 0$ of a system in thermal equilibrium with a reservoir the density matrix get close to the ground state:
$$
\rho=\frac{e^{-\beta H}}{tr(e^{-\beta H})}\,\,\,\,\rightarrow \,\,\,\,|0\rangle\langle 0|
$$
And then you can use the path integral ideas seeing the $\beta$ as a time, and the limit as a $\beta\rightarrow \infty$.
