  Many approximate solutions of the time-dependent Schr\"odinger equation can
be formulated as exact solutions of a nonlinear Schr\"odinger equation with an
effective Hamiltonian operator depending on the state of the system. We show
that several well-known Gaussian wavepacket dynamics methods, such as Heller's
original thawed Gaussian approximation or Coalson and Karplus's variational
Gaussian approximation, fit into this framework if the effective potential is a
quadratic polynomial with state-dependent coefficients. We study such a
nonlinear Schr\"odinger equation in general: in particular, we derive general
equations of motion for the Gaussian's parameters, demonstrate the time
reversibility and norm conservation, and analyze conservation of the energy,
effective energy, and symplectic structure. We also describe efficient
geometric integrators of arbitrary even orders of accuracy in the time step for
the numerical solution of this nonlinear Schr\"odinger equation. The general
presentation is illustrated by examples of this family of Gaussian wavepacket
dynamics, including the variational and nonvariational thawed and frozen
Gaussian approximations, and special limits of these methods based on the
global harmonic, local harmonic, single-Hessian, local cubic, and single
quartic approximations for the potential energy. Without substantially
increasing the cost, the proposed single quartic variational thawed Gaussian
wavepacket dynamics improves the accuracy over the local cubic approximation
and, at the same time, conserves both the effective energy and symplectic
structure, in contrast to the much more expensive local quartic approximation.
Most results are presented in both Heller's and Hagedorn's parametrizations of
the Gaussian wavepacket.
