Not yet reviewed by Pith; the record is open.

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.