  Geometric integrators of the Schr\"{o}dinger equation conserve exactly many
invariants of the exact solution. Among these integrators, the split-operator
algorithm is explicit and easy to implement, but, unfortunately, is restricted
to systems whose Hamiltonian is separable into a kinetic and potential terms.
Here, we describe several implicit geometric integrators applicable to both
separable and non-separable Hamiltonians, and, in particular, to the
nonadiabatic molecular Hamiltonian in the adiabatic representation. These
integrators combine the dynamic Fourier method with recursive symmetric
composition of the trapezoidal rule or implicit midpoint method, which results
in an arbitrary order of accuracy in the time step. Moreover, these integrators
are exactly unitary, symplectic, symmetric, time-reversible, and stable, and,
in contrast to the split-operator algorithm, conserve energy exactly,
regardless of the accuracy of the solution. The order of convergence and
conservation of geometric properties are proven analytically and demonstrated
numerically on a two-surface NaI model in the adiabatic representation.
Although each step of the higher order integrators is more costly, these
algorithms become the most efficient ones if higher accuracy is desired; a
thousand-fold speedup compared to the second-order trapezoidal rule (the
Crank-Nicolson method) was observed for wavefunction convergence error of
$10^{-10}$. In a companion paper [J. Roulet, S. Choi, and J. Van\'{\i}\v{c}ek
(2019)], we discuss analogous, arbitrary-order compositions of the
split-operator algorithm and apply both types of geometric integrators to a
higher-dimensional system in the diabatic representation.
