  The viability of a variant of numerical stochastic perturbation theory, where
the Langevin equation is replaced by the SMD algorithm, is examined. In
particular, the convergence of the process to a unique stationary state is
rigorously established and the use of higher-order symplectic integration
schemes is shown to be highly profitable in this context. For illustration, the
gradient-flow coupling in finite volume with Schr\"odinger functional boundary
conditions is computed to two-loop (i.e. NNL) order in the SU(3) gauge theory.
The scaling behaviour of the algorithm turns out to be rather favourable in
this case, which allows the computations to be driven close to the continuum
limit.
