  We present a new method for sampling stochastic displacements in Brownian
Dynamics (BD) simulations of colloidal scale particles. The method relies on a
new formulation for Ewald summation of the Rotne-Prager-Yamakawa (RPY) tensor,
which guarantees that the real-space and wave-space contributions to the tensor
are independently symmetric and positive-definite for all possible particle
configurations. Brownian displacements are drawn from a superposition of two
independent samples: a wave-space (far-field or long-ranged) contribution,
computed using techniques from fluctuating hydrodynamics and non-uniform Fast
Fourier Transforms; and a real-space (near-field or short-ranged) correction,
computed using a Krylov subspace method. The combined computational complexity
of drawing these two independent samples scales linearly with the number of
particles. The proposed method circumvents the super-linear scaling exhibited
by all known iterative sampling methods applied directly to the RPY tensor that
results from the power law growth of the condition number of tensor with the
number of particles. Calculations for hard sphere dispersions and colloidal
gels are illustrated and used to explore the role of microstructure on
performance of the algorithm. In practice, the logarithmic part of the
predicted scaling is not observed and the algorithm scales linearly for up to 4
million particles, obtaining speed ups of over an order of magnitude over
existing iterative methods, and making the cost of computing Brownian
displacements comparable the cost of computing deterministic displacements in
BD simulations. A high-performance implementation employing non-uniform fast
Fourier transforms implemented on graphics processing units and integrated with
the software package HOOMD-blue is used for benchmarking.
