  In light of the recently developed complete GJ set of single random variable
stochastic, discrete-time St{\o}rmer-Verlet algorithms for statistically
accurate simulations of Langevin equations, we investigate two outstanding
questions: 1) Are there any algorithmic or statistical benefits from including
multiple random variables per time-step, and 2) are there objective reasons for
using one or more methods from the available set of statistically correct
algorithms? To address the first question, we assume a general form for the
discrete-time equations with two random variables and then follow the
systematic, brute-force GJ methodology by enforcing correct thermodynamics in
linear systems. It is concluded that correct configurational Boltzmann sampling
of a particle in a harmonic potential implies correct configurational
free-particle diffusion, and that these requirements only can be accomplished
if the two random variables per time step are identical. We consequently submit
that the GJ set represents all possible stochastic St{\o}rmer-Verlet methods
that can reproduce time-step-independent statistics of linear systems. The
second question is thus addressed within the GJ set. Based in part on numerical
simulations of complex molecular systems, and in part on analytic scaling of
time, we analyze the apparent difference in stability between different
methods. We attribute this difference to the inherent time scaling in each
method, and suggest that this scaling may lead to inconsistencies in the
interpretation of dynamical and statistical simulation results. We therefore
suggest that the method with the least inherent time-scaling, the GJ-I/GJF-2GJ
method, be preferred for statistical applications where spurious rescaling of
time is undesirable.
