^{1}, Edmond Chow

^{2}, Yousef Saad

^{3}and Jeffrey Skolnick

^{1}

### 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.

This research was supported in part by National of Institutes of Health (NIH) Grant No. GM-37408 of the Division of General Medical Sciences of the National Institutes of Health. Dr. Saad was supported by National Science Foundation (NSF) Grant No. NSF/DMS-0810938.

I. INTRODUCTION

II. THEORY

A. Chebyshev polynomial approximations

B. Truncated expansion approximation

C. Krylov subspace approximation

1. Lanczos process

2. Block-Lanczos process for multiple vectors

III. MODELS AND SIMULATION METHODS

A. Random polymer chain model

B. Collapsed chain model

C. Monodisperse suspension model

D. Brownian dynamics simulation and analysis

IV. RESULTS AND DISCUSSION

A. Convergence

1. Convergence rate and error estimates

2. Effect of block size on convergence

3. Simulation model properties affecting convergence

B. BD simulations

1. Update interval of diffusiontensor

2. Required accuracy of Brownian noise vectors

3. Comparison of covariance matrix generated from Brownian noise vectors

C. Computational time

V. CONCLUSIONS

### Key Topics

- Subspaces
- 79.0
- Polymers
- 48.0
- Eigenvalues
- 32.0
- Polynomials
- 32.0
- Tensor methods
- 17.0

##### C08

## Figures

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.

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.

(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.

(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.

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.

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.

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.

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.

Values of *E* _{1} 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.

Values of *E* _{1} 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.

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.

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.

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.

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.

## Tables

Krylov subspace algorithm for computing

Krylov subspace algorithm for computing

Block-Krylov subspace algorithm for computing a block of *s* correlated vectors, , each vector with distribution *N*(0, **D**).

Block-Krylov subspace algorithm for computing a block of *s* correlated vectors, , each vector with distribution *N*(0, **D**).

Errors (%) in *D* obtained with various diffusion matrix update intervals, *λ* _{RPY}, relative to results with *λ* _{RPY} = 1 for the monodisperse model at various volume fractions Φ with *N* = 200.

Errors (%) in *D* obtained with various diffusion matrix update intervals, *λ* _{RPY}, relative to results with *λ* _{RPY} = 1 for the monodisperse model at various volume fractions Φ with *N* = 200.

Errors (%) in *D* _{cm} obtained from the simulations using the Krylov method with various levels of Brownian noise accuracies and the TEA method relative to results using the Cholesky method with *λ* _{RPY} = 1 for various chain lengths for the random polymer model.

Errors (%) in *D* _{cm} obtained from the simulations using the Krylov method with various levels of Brownian noise accuracies and the TEA method relative to results using the Cholesky method with *λ* _{RPY} = 1 for various chain lengths for the random polymer model.

Errors (%) in *D* _{cm} at equilibrated states obtained from the simulations using the Krylov method with various Brownian noise accuracies and the TEA method relative to results using the Cholesky method with *λ* _{RPY} = 1 for various chain lengths for the collapsed polymer model.

Errors (%) in *D* _{cm} at equilibrated states obtained from the simulations using the Krylov method with various Brownian noise accuracies and the TEA method relative to results using the Cholesky method with *λ* _{RPY} = 1 for various chain lengths for the collapsed polymer model.

Errors (%) in *D* obtained from the simulations using the Krylov method with various Brownian noise accuracies and the TEA method relative to results using the Cholesky method with *λ* _{RPY} = 1 for monodisperse model at various volume fractions Φ with number of particles *N* = 200.

Errors (%) in *D* obtained from the simulations using the Krylov method with various Brownian noise accuracies and the TEA method relative to results using the Cholesky method with *λ* _{RPY} = 1 for monodisperse model at various volume fractions Φ with number of particles *N* = 200.

Article metrics loading...

Full text loading...

Commenting has been disabled for this content