Many approximate solutions of the time-dependent Schrodinger equation can be formulated as exact solutions of a nonlinear Schrodinger equation with an effective Hamiltonian operator depending on the state of the system. We show that Heller's thawed Gaussian approximation, Coalson and Karplus's variational Gaussian approximation, and other Gaussian wavepacket dynamics methods fit into this framework if the effective potential is a quadratic polynomial with state-dependent coefficients. We study such a nonlinear Schrodinger equation in full generality: we derive general equations of motion for the Gaussian's parameters, demonstrate time reversibility and norm conservation, and analyze conservation of energy, effective energy, and symplectic structure. We also describe efficient, high-order geometric integrators for the numerical solution of this nonlinear Schrodinger equation. The general theory is illustrated by examples of this family of Gaussian wavepacket dynamics, including the variational and nonvariational thawed and frozen Gaussian approximations and their special limits based on the global harmonic, local harmonic, single-Hessian, local cubic, and local quartic approximations for the potential energy. We also propose a new method by augmenting the local cubic approximation with a single fourth derivative. Without substantially increasing the cost, the proposed "single-quartic" variational Gaussian approximation improves the accuracy over the local cubic approximation and, at the same time, conserves both the effective energy and symplectic structure, unlike the much more expensive local quartic approximation. Most results are presented in both Heller's and Hagedorn's parametrizations of the Gaussian wavepacket.