The nuclear-electronic orbital (NEO) method for the calculation of mixed nuclear-electronic wave functions is presented. Both electronic and nuclear molecular orbitals are expressed as linear combinations of Gaussian basis functions. In the NEO-HF (Hartree-Fock) method, the energy corresponding to the single-configurational mixed nuclear-electronic wave function is minimized with respect to the molecular orbitals. Multiconfigurational approaches are implemented to include significant correlation effects. In the NEO-CI (configuration interaction) method, the energy corresponding to the multiconfigurational mixed nuclear-electronic wave function is minimized with respect to the CI coefficients. In the NEO-MCSCF (multiconfigurational self-consistent-field) method, the energy is minimized with respect to the molecular orbitals as well as the CI coefficients. Analytic gradient expressions are presented for NEO-HF and NEO-MCSCF. These analytic gradients allow the variational optimization of the centers of the nuclear basis functions. They also enable the location and characterization of geometry stationary points and the generation of minimum energy paths and dynamic reaction paths. The advantages of the NEO approach are that nuclear quantum effects are incorporated during the electronic structure calculation, the Born-Oppenheimer separation of electrons and nuclei is avoided, excited vibrational-electronic states may be calculated, and its accuracy may be improved systematically. Initial applications are presented to illustrate the computational feasibility and accuracy of this approach.