The multiconfiguration time-dependent Hartree-Fock (MCTDHF) method is formulated for treating the coupled electronic and nuclear dynamics of diatomic molecules without the Born- Oppenheimer approximation. The method treats the full dimensionality of the electronic motion, uses no model interactions, and is in principle capable of an exact nonrelativistic description of diatomics in electromagnetic fields. An expansion of the wave function in terms of configurations of orbitals whose dependence on internuclear distance is only that provided by the underlying prolate spheroidal coordinate system is demonstrated to provide the key simplifications of the working equations that allow their practical solution. Photoionization cross sections are also computed from the MCTDHF wave function in calculations using short pulses.