  Application of the time-dependent variational principle to a linear
combination of frozen-width Gaussians describing the nuclear wavefunction
provides a formalism where the total energy is conserved. The computational
downside of this formalism is that trajectories of individual Gaussians are
solutions of a coupled system of differential equations, limiting
implementation to serial propagation algorithms. To allow for parallelization
and acceleration of the computation, independent trajectories based on
simplified equations of motion were suggested. Unfortunately, within practical
realizations involving finite Gaussian bases, this simplification leads to
breaking the energy conservation. We offer a solution for this problem by using
Lagrange multipliers to ensure the energy and norm conservation regardless of
basis function trajectories or basis completeness. We illustrate our approach
within the Multi-Configuration Ehrenfest method considering a linear vibronic
coupling model.
