  Exact nonadiabatic quantum evolution preserves many geometric properties of
the molecular Hilbert space. In a companion paper [S. Choi and J.
Van\'{\i}\v{c}ek, 2019], we presented numerical integrators of arbitrary-order
of accuracy that preserve these geometric properties exactly even in the
adiabatic representation, in which the molecular Hamiltonian is not separable
into a kinetic and potential terms. Here, we focus on the separable Hamiltonian
in diabatic representation, where the split-operator algorithm provides a
popular alternative because it is explicit and easy to implement, while
preserving most geometric invariants. Whereas the standard version has only
second-order accuracy, we implemented, in an automated fashion, its recursive
symmetric compositions, using the same schemes as in the companion paper, and
obtained integrators of arbitrary even order that still preserve the geometric
properties exactly. Because the automatically generated splitting coefficients
are redundant, we reduce the computational cost by pruning these coefficients
and lower memory requirements by identifying unique coefficients. The order of
convergence and preservation of geometric properties are justified analytically
and confirmed numerically on a one-dimensional two-surface model of NaI and a
three-dimensional three-surface model of pyrazine. As for efficiency, we find
that to reach a convergence error of 10$^{-10}$, a 600-fold speedup in the case
of NaI and a 900-fold speedup in the case of pyrazine are obtained with the
higher-order compositions instead of the second-order split-operator algorithm.
The pyrazine results suggest that the efficiency gain survives in higher
dimensions.
