Fast algorithms for Brownian dynamics simulation with hydrodynamic interactions

In the Brownian dynamics simulation with hydrodynamic interactions, one needs to generate the total displacement vectors of Brownian particles consisting of two parts: a deterministic part which is proportional to the product of the Rotne-Prager-Yamakawa (RPY) tensor D [1, 2] and the given external forces F; and a hydrodynamically correlated random part whose covariance is proportional to the RPY tensor. To be more precise, one needs to calculate Du for a given vector u and compute √Dv for a normally distributed random vector v. For an arbitrary N-particle configuration, D is a 3N x 3N matrix and u, v are vectors of length 3N. Thus, classical algorithms require O(N2) operations for computing Du and O(N3) operations for computing √Dv, which are prohibitively expensive and render large scale simulations impossible since one needs to carry out these calculations many times in a Brownian dynamics simulation. In this dissertation, we first present two fast multipole methods (FMM) for computing Du. The first FMM is a simple application of the kernel independent FMM (KIFMM) developed by Ying, Biros, and Zorin [3], which requires 9 scalar FMM calls. The second FMM, similar to the FMM for Stokeslet developed by Tornberg and Greengard [4], decomposes the RPY tensor into harmonic potentials and its derivatives, and thus requires only four harmonic FMM calls. Both FMMs reduce the computational cost of Du from O(N2) to O(N) for an arbitrary N-particle configuration. We then discuss several methods of computing √Dv, which are all based on the Krylov subspace approximations, that is, replacing √Dv by p(D)v with p(D) a low degree polynomial in D. We first show rigorously that the popular Chebyshev spectral approximation method (see, for example, [5, 6]) requires √κ log 1/ε terms for a desired precision E, where K is the condition number of the RPY tensor D. In the Chebyshev spectral approximation method, one also needs to estimate the extreme eigenvalues of D. We have considered several methods: the classical Lanczos method, the Chebyshev-Davidson method, and the safeguarded Lanczos method proposed by Zhou and Li [7]. Our numerical experiments indicate that K is usually very small when the particles are distributed uniformly with low density, and that the safeguarded Lanczos method is most effective for our cases with very little additional computational cost. Thus, when combined with the FMMs we described earlier, the Chebyshev approximation method with safeguarded Lanczos method as eigenvalue estimators essentially reduces the cost of computing √Dv from O(N3) to O(N) for most practical particle configurations. Finally, we propose to combine the so-called spectral Lanczos decomposition method (SLDM) (see, for example, [8]) and the FMMs to compute √Dv. Our numerical experiments show that the SLDM is generally more efficient than the popular Chebyshev spectral approximation method. The fast algorithms developed in this dissertation will be useful for the study of diffusion limited reactions, polymer dynamics, protein folding, and particle coagulation as it enables large scale Brownian dynamics simulations. Moreover, the algorithms can be extended to speed up the computation involving the matrix square root for many other matrices, which has potential applications in areas such as statistical analysis with certain spatial correlations and model reduction in dynamic control theory

Hydrodynamics of Suspensions of Passive and Active Rigid Particles: A Rigid Multiblob Approach

We develop a rigid multiblob method for numerically solving the mobility problem for suspensions of passive and active rigid particles of complex shape in Stokes flow in unconfined, partially confined, and fully confined geometries. As in a number of existing methods, we discretize rigid bodies using a collection of minimally-resolved spherical blobs constrained to move as a rigid body, to arrive at a potentially large linear system of equations for the unknown Lagrange multipliers and rigid-body motions. Here we develop a block-diagonal preconditioner for this linear system and show that a standard Krylov solver converges in a modest number of iterations that is essentially independent of the number of particles. For unbounded suspensions and suspensions sedimented against a single no-slip boundary, we rely on existing analytical expressions for the Rotne-Prager tensor combined with a fast multipole method or a direct summation on a Graphical Processing Unit to obtain an simple yet efficient and scalable implementation. For fully confined domains, such as periodic suspensions or suspensions confined in slit and square channels, we extend a recently-developed rigid-body immersed boundary method to suspensions of freely-moving passive or active rigid particles at zero Reynolds number. We demonstrate that the iterative solver for the coupled fluid and rigid body equations converges in a bounded number of iterations regardless of the system size. We optimize a number of parameters in the iterative solvers and apply our method to a variety of benchmark problems to carefully assess the accuracy of the rigid multiblob approach as a function of the resolution. We also model the dynamics of colloidal particles studied in recent experiments, such as passive boomerangs in a slit channel, as well as a pair of non-Brownian active nanorods sedimented against a wall.Comment: Under revision in CAMCOS, Nov 201

Fast Algorithms for Brownian Dynamics with Hydrodynamic Interactions

In this dissertation, we contributed on three fundamental parts of Brownian dynamics simulations with hydrodynamic interactions. The first part of the dissertation is to derive the formulas for computing the electric field gradients by the new version of fast multipole method(FMM) and to implement them as new functions for existing FMM solvers. In the second part of the dissertation, we discuss how to decompose the far-field Rotne-Prager-Yamakawa potential into four far-field Laplace FMM calls including electrostatic potential, electric field and field gradient terms. A parallelized Rotne-Prager-Yamakawa solver based on the new version of fast multipole method has been developed with tunable accuracy. The solver makes it computationally viable for large-scale, long-time Brownian dynamic simulations with hydrodynamic interactions. In the third part, a model is built toward an accurate description of hydrodynamic effects on the translational and rotational dynamics of complex, rigid macromolecules with arbitrary shape in suspension. The grand diffusion matrix is calculated by employing the bead-shell model for describing the shape and structure of macromolecules in the many-body system. Two fast algorithms based on block conjugate gradient method and the Schur complement method are developed for computing the translational and angular velocities, as well as the displacements and orientations in order to track the trajectories of the macromolecules in the complex structured biological system.Doctor of Philosoph

A general formulation of Bead Models applied to flexible fibers and active filaments at low Reynolds number

This contribution provides a general framework to use Lagrange multipliers for the simulation of low Reynolds number fiber dynamics based on Bead Models (BM). This formalism provides an efficient method to account for kinematic constraints. We illustrate, with several examples, to which extent the proposed formulation offers a flexible and versatile framework for the quantitative modeling of flexible fibers deformation and rotation in shear flow, the dynamics of actuated filaments and the propulsion of active swimmers. Furthermore, a new contact model called Gears Model is proposed and successfully tested. It avoids the use of numerical artifices such as repulsive forces between adjacent beads, a source of numerical difficulties in the temporal integration of previous Bead Models.Comment: 41 pages, 15 figure

Hierarchical Orthogonal Matrix Generation and Matrix-Vector Multiplications in Rigid Body Simulations

In this paper, we apply the hierarchical modeling technique and study some numerical linear algebra problems arising from the Brownian dynamics simulations of biomolecular systems where molecules are modeled as ensembles of rigid bodies. Given a rigid body $p$ consisting of $n$ beads, the $6 \times 3n$ transformation matrix $Z$ that maps the force on each bead to $p$'s translational and rotational forces (a $6\times 1$ vector), and $V$ the row space of $Z$, we show how to explicitly construct the $(3n-6) \times 3n$ matrix $\tilde{Q}$ consisting of $(3n-6)$ orthonormal basis vectors of $V^{\perp}$ (orthogonal complement of $V$) using only $\mathcal{O}(n \log n)$ operations and storage. For applications where only the matrix-vector multiplications $\tilde{Q}{\bf v}$ and $\tilde{Q}^T {\bf v}$ are needed, we introduce asymptotically optimal $\mathcal{O}(n)$ hierarchical algorithms without explicitly forming $\tilde{Q}$. Preliminary numerical results are presented to demonstrate the performance and accuracy of the numerical algorithms