Krylov subspace methods for computing hydrodynamic interactions in brownian dynamics simulations

Abstract

Hydrodynamic interactions play an important role in the dynamics of macromolecules. The most common way to take into account hydrodynamic effects in molecular simulations is in the context of a brownian dynamics simulation. However, the calculation of correlated brownian noise vectors in these simulations is computationally very demanding and alternative methods are desirable. This paper studies methods based on Krylov subspaces for computing brownian noise vectors. These methods are related to Chebyshev polynomial approximations, but do not require eigenvalue estimates. We show that only low accuracy is required in the brownian noise vectors to accurately compute values of dynamic and static properties of polymer and monodisperse suspension models. With this level of accuracy, the computational time of Krylov subspace methods scales very nearly as O(N(2)) for the number of particles N up to 10 000, which was the limit tested. The performance of the Krylov subspace methods, especially the "block" version, is slightly better than that of the Chebyshev method, even without taking into account the additional cost of eigenvalue estimates required by the latter. Furthermore, at N = 10,000, the Krylov subspace method is 13 times faster than the exact Cholesky method. Thus, Krylov subspace methods are recommended for performing large-scale brownian dynamics simulations with hydrodynamic interactions.

Figure 1

Convergence of (a) the Krylov subspace method and (b) the Chebyshev method measured by various error estimates, E_k^exact, E_k, and E_f. Diffusion matrices were constructed from five configurations of random polymer chains with length N = 1000 and the results represent the average over the five configurations. Standard deviations for all data points are so small that they are not displayed. E_k with k = 1 was set to 1.

Figure 2

(a) Effect of block size on convergence rate in the block Krylov subspace method. The errors are computed for the first vector of the block of vectors. Results are averages of the five random polymer chains with N = 1000. (b) Comparison of E_k^exact and E_k in the block Krylov subspace method. Results for polymers with N = 1000 and block size = 50 are shown. All results are the average over the five independent configurations. Standard deviations for all data points are so small that they are not displayed. E_k with k = 1 was set to 1.

Figure 3

Effect of model, number of particles, and volume fraction on convergence of the Krylov subspace method. (a) Convergence of the random and collapsed polymer models with N = 200 and 1000. (b) Convergence of the monodisperse suspension model at various volume fractions Φ = 0.1, 0.2, 0.3, 0.4, and 0.5 with N = 200 and N = 1000. For monodisperse suspensions, convergence is slower for larger volume fractions and for larger numbers of particles.

Figure 4

Effect of λ_RPY on dynamic properties for the random and collapsed polymers, and monodisperse suspension model. (a) D_cm for the random and collapsed polymers with various polymer lengths N obtained from BD simulation with various λ_RPY. Lines are fit to the data for λ_RPY = 1 with D_cm ∝ N^−ν and their exponents are shown. For both polymer models, results with different λ_RPY are so close that their plots overlap and are almost indistinguishable in the figure. (b) D for the monodisperse suspension model with number of particles of 200 at various volume fractions Φ obtained from BD simulations with various λ_RPY. The results with λ_RPY = 1 are connected by a broken line to guide the eye. The Cholesky factorization method was used in the BD simulations.

Figure 5

Values of E₁ for the covariance matrices constructed from sets of Brownian noise vectors generated by the Cholesky, Krylov subspace with E_k = 0.1, and TEA methods. Results of the random polymer chain model (a, d), the collapsed polymer model (b, e), and the monodisperse suspension model at volume fraction of 0.3 (c, f) are shown. Left subfigures (a, b, c) are for N = 200 and right subfigures (d, e, f) are for N = 1000. For the monodisperse suspension model, results of the Cholesky (red lines) and Krylov subspace with E_k = 0.1 (green lines) are so close that their lines overlap and are indistinguishable in the figure. All results are the average over five independent configurations. Standard deviations for all data points are so small that they are not displayed.

Figure 6

Computational time required for generating 100 correlated Brownian noise vectors by the block-Krylov subspace method with different block sizes. Timings for (a) E_k = 0.1 and (b) E_k = 0.01 are shown. The random polymer model with length N = 1000–10 000 was used for this timing test. For block size = 1, the DSYMV BLAS routine was employed for matrix-vector multiplications. For other block sizes, the DSYMM BLAS routine was employed for matrix-matrix multiplications. The latter routine is not optimized for matrix-vector multiplications.

Figure 7

Scaling of number of iterations and computational time of the block-Krylov subspace method with the number of particles. The random polymer model with length N = 1000–10 000 was used for this timing test. Number of iterations required with thresholds (a) E_k = 0.1 and (b) E_k = 0.01 for the block-Krylov subspace and Cholesky methods. Computational time required for generating 50 correlated Brownian noise vectors by the block-Krylov subspace and Cholesky methods with (c) E_k = 0.1 and (d) E_k = 0.01. A block size of 50 was used for the block-Krylov subspace method. Results for the Cholesky and TEA methods are also shown for comparison (also computed in block fashion). Dashed lines are fitted linear slopes for a range of N = 4000–10 000. The values of slopes for these methods are shown in inside the figures.