  We study the energy levels of H$_2$ molecules in a superstrong magnetic field
($B\go 10^{12}$ G), typically found on the surfaces of neutron stars. The
interatomic interaction potentials are calculated by a Hartree-Fock method with
multi-configurations assuming electrons are in the ground Landau state. Both
the aligned configurations and arbitrary orientations of the molecular axis
with respect to the magnetic field axis are considered. Different types of
molecular excitations are then studied: electronic excitations, aligned (along
the magnetic axis) vibrational excitations, transverse vibrational excitations
(a constrained rotation of the molecular axis around the magnetic field line).
Similar results for the molecular ion H$_2^+$ are also obtained and compared
with previous variational calculations. Both numerical results and analytical
fitting formulae are given for a wide range of field strengths. In contrast to
the zero-field case, it is found that the transverse vibrational excitation
energies can be larger than the aligned vibration excitation, and they both can
be comparable or larger than the electronic excitations. For $B\go
B_{crit}=4.23\times 10^{13}$ G, the Landau energy of proton is appreciable and
there is some controversy regarding the dissociation energy of H$_2$. We show
that H$_2$ is bound even for $B>>B_{crit}$ and that neither proton has a Landau
excitation in the ground molecular state.
