An immersed boundary method for Brownian dynamics simulation of polymers in complex geometries: application to DNA flowing through a nanoslit with embedded nanopits

This work presents an immersed boundary method that allows fast Brownian dynamics simulation of solutions of polymer chains and other Brownian objects in complex geometries with fluctuating hydrodynamics. The approach is based on the general geometry Ewald-like method, which solves the Stokes equation with distributed regularized point forces in O(N) or O(NlogN) operations, where N is the number of point forces in the system. Time-integration is performed using a midpoint algorithm and Chebyshev polynomial approximation proposed by Fixman. This approach is applied to the dynamics of a genomic DNA molecule driven by flow through a nanofluidic slit with an array of nanopits on one wall of the slit. The dynamics of the DNA molecule was studied as a function of the Péclet number and chain length (the base case being λ-DNA). The transport characteristics of the hopping dynamics in this device differ at low and high Péclet number, and for long DNA, relative to the pit size, the dynamics is governed by the segments residing in the pit. By comparing with results that neglect them, hydrodynamic interactions are shown to play an important quantitative role in the hopping dynamics.     

Pair diffusion, hydrodynamic interactions, and available volume in dense fluids

We calculate the pair diffusion coefficient D(r) as a function of the distance r between two hard sphere particles in a dense monodisperse fluid. The distance-dependent pair diffusion coefficient describes the hydrodynamic interactions between particles in a fluid that are central to theories of polymer and colloid dynamics. We determine D(r) from the propagators (Green's functions) of particle pairs obtained from molecular dynamics simulations. At distances exceeding ~3 molecular diameters, the calculated pair diffusion coefficients are in excellent agreement with predictions from exact macroscopic hydrodynamic theory for large Brownian particles suspended in a solvent bath, as well as the Oseen approximation. However, the asymptotic 1/r distance dependence of D(r) associated with hydrodynamic effects emerges only after the pair distance dynamics has been followed for relatively long times, indicating non-negligible memory effects in the pair diffusion at short times. Deviations of the calculated D(r) from the hydrodynamic models at short distances r reflect the underlying many-body fluid structure, and are found to be correlated to differences in the local available volume. The procedure used here to determine the pair diffusion coefficients can also be used for single-particle diffusion in confinement with spherical symmetry.     

Krylov subspace methods for computing hydrodynamic interactions in Brownian dynamics simulations

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.     

Brownian dynamics simulations of a flexible polymer chain which includes continuous resistance and multibody hydrodynamic interactions

Using methods adapted from the simulation of suspension dynamics, we have developed a Brownian dynamics algorithm with multibody hydrodynamic interactions for simulating the dynamics of polymer molecules. The polymer molecule is modeled as a chain composed of a series of inextensible, rigid rods with constraints at each joint to ensure continuity of the chain. The linear and rotational velocities of each segment of the polymer chain are described by the slender-body theory of Batchelor [J. Fluid Mech. 44, 419 (1970)]. To include hydrodynamic interactions between the segments of the chain, the line distribution of forces on each segment is approximated by making a Legendre polynomial expansion of the disturbance velocity on the segment, where the first two terms of the expansion are retained in the calculation. Thus, the resulting linear force distribution is specified by a center of mass force, couple, and stresslet on each segment. This method for calculating the hydrodynamic interactions has been successfully used to simulate the dynamics of noncolloidal suspensions of rigid fibers [O. G. Harlen, R. R. Sundararajakumar, and D. L. Koch, J. Fluid Mech. 388, 355 (1999); J. E. Butler and E. S. G. Shaqfeh, J. Fluid Mech. 468, 204 (2002)]. The longest relaxation time and center of mass diffusivity are among the quantities calculated with the simulation technique. Comparisons are made for different levels of approximation of the hydrodynamic interactions, including multibody interactions, two-body interactions, and the "freely draining" case with no interactions. For the short polymer chains studied in this paper, the results indicate a difference in the apparent scaling of diffusivity with polymer length for the multibody versus two-body level of approximation for the hydrodynamic interactions.     

Diffusion and velocity relaxation of a Brownian particle immersed in a viscous compressible fluid confined between two parallel plane walls

The diffusion tensor and velocity correlation function of a Brownian particle immersed in a viscous compressible fluid confined between two parallel plane walls are calculated in point approximation. The fluid is assumed to satisfy stick boundary conditions at the walls. It is found that the velocity correlation function decays asymptotically with a negative t(-2) long-time tail due to coupling to overdamped sound waves. The coefficient of the long-time tail is calculated and shown to be independent of fluid viscosity.     

Brownian dynamics simulations of the self- and collective rotational diffusion coefficients of rigid long thin rods

Recently a microscopic theory for the dynamics of suspensions of long thin rigid rods was presented, confirming and expanding the well-known theory by Doi and Edwards [The Theory of Polymer Dynamics (Clarendon, Oxford, 1986)] and Kuzuu [J. Phys. Soc. Jpn. 52, 3486 (1983)]. Here this theory is put to the test by comparing it against computer simulations. A Brownian dynamics simulation program was developed to follow the dynamics of the rods, with a length over a diameter ratio of 60, on the Smoluchowski time scale. The model accounts for excluded volume interactions between rods, but neglects hydrodynamic interactions. The self-rotational diffusion coefficients D(r)(phi) of the rods were calculated by standard methods and by a new, more efficient method based on calculating average restoring torques. Collective decay of orientational order was calculated by means of equilibrium and nonequilibrium simulations. Our results show that, for the currently accessible volume fractions, the decay times in both cases are virtually identical. Moreover, the observed decay of diffusion coefficients with volume fraction is much quicker than predicted by the theory, which is attributed to an oversimplification of dynamic correlations in the theory.     

Equilibrium conformational dynamics of a polymer in a solvent

Molecular dynamics simulations were used to study the conformational dynamics of a bead-spring model polymer in an explicit solvent under good solvent conditions. The dynamics of the polymer chain were investigated using an analysis of the time autocorrelation functions of the Rouse coordinates of the polymer chain. We have investigated the variation of the correlation functions with polymer chain length N, solvent density rho, and system size. The measured initial decay rates gamma(p) of the correlation functions were compared with the predictions from a theory of polymer dynamics which uses the Oseen tensor to describe hydrodynamic interactions between monomers. Over the range of chain lengths considered (N = 30-60 monomers), the predicted scaling of gamma(p) proportional to N(-3nu) was observed at high rho, where nu is the polymer scaling exponent. The predicted gamma(p) are generally higher than the measured values. This discrepancy increases with decreasing rho, as a result in the breakdown in the conditions required for the Oseen approximation. The agreement between theory and simulation at high rho improves considerably if the theoretical expression for gamma(p) is modified to avoid sum-to-integral approximations, and if the values of (R(p)2), which are used in the theory, are taken directly from the simulation rather than being calculated using approximate scaling relations. The observed finite-size scaling of gamma(p) is not quantitatively consistent with the theoretical predictions.     

Polymer chains in a soft nanotube: a Monte Carlo study

We study the equilibrium properties of flexible polymer chains confined in a soft tube by means of extensive Monte Carlo simulations. The tube wall is that of a single sheet six-coordinated self-avoiding tethered membrane. Our study assumes that there is no adsorption of the chain on the wall. By varying the length N of the polymer and the tube diameter D we examine the variation of the polymer gyration radius Rg and diffusion coefficient Ddiff in soft and rigid tubes of identical diameter and compare them to scaling theory predictions. We find that the swollen region of the soft tube surrounding the chain exhibits a cigarlike cylindrical shape for sufficiently narrow tubes with D相似文献   

Conformation and diffusion behavior of ring polymers in solution: a comparison between molecular dynamics, multiparticle collision dynamics, and lattice Boltzmann simulations

We have studied the effect of chain topology on the structural properties and diffusion of polymers in a dilute solution in a good solvent. Specifically, we have used three different simulation techniques to compare the chain size and diffusion coefficient of linear and ring polymers in solution. The polymer chain is modeled using a bead-spring representation. The solvent is modeled using three different techniques: molecular dynamics (MD) simulations with a particulate solvent in which hydrodynamic interactions are accounted through the intermolecular interactions, multiparticle collision dynamics (MPCD) with a point particle solvent which has stochastic interactions with the polymer, and the lattice Boltzmann method in which the polymer chains are coupled to the lattice fluid through friction. Our results show that the three methods give quantitatively similar results for the effect of chain topology on the conformation and diffusion behavior of the polymer chain in a good solvent. The ratio of diffusivities of ring and linear polymers is observed to be close to that predicted by perturbation calculations based on the Kirkwood hydrodynamic theory.     

Strong effect of hydrodynamic coupling on the electric dichroism of bent rods

The effect of hydrodynamic coupling on the spatial orientation of rigid bent rods in electric fields has been analyzed by Brownian dynamics simulations. Bead models for smoothly bent rods were constructed with dimensions of DNA double helices, and established simulation procedures were used to calculate their diffusion tensor, including the translational-rotational coupling tensor. The electric and optical parameters were assigned on the basis of known properties of double helices. Brownian dynamics simulations of the orientation of these models in electric fields showed that both transients and amplitudes of the calculated dichroism are very strongly dependent on translational-rotational coupling over a wide range of electric field strengths. For example, the stationary dichroism of a smoothly bent 179 bp DNA fragment calculated at low field strengths is positive in the presence and negative in the absence of hydrodynamic coupling. The transients are converted from a biphasic to a monophasic shape, when hydrodynamic coupling is turned off. The large changes resulting from hydrodynamic coupling were controlled by calculations based on analytical expressions derived for electrooptical response curves in the limit of low electric field strengths; the results obtained by this independent approach are in very satisfactory agreement with our Brownian dynamics simulations. The effect is strongly dependent on the electric dipole and on its direction. In the absence of any dipole the coupling effect was not observed. The coupling effect increases with the size of the bent rods. Because most macromolecular structures are known to have induced and/or permanent dipole moments, large effects of hydrodynamic coupling on both the amplitudes and the transients of the electric dichroism/birefringence must be expected in general for structures with nonsymmetric shape.     

Lattice-Boltzmann simulations of the dynamics of polymer solutions in periodic and confined geometries

A numerical method to simulate the dynamics of polymer solutions in confined geometries has been implemented and tested. The method combines a fluctuating lattice-Boltzmann model of the solvent [Ladd, Phys. Rev. Lett. 70, 1339 (1993)] with a point-particle model of the polymer chains. A friction term couples the monomers to the fluid [Ahlrichs and Dunweg, J. Chem. Phys. 111, 8225 (1999)], providing both the hydrodynamic interactions between the monomers and the correlated random forces. The coupled equations for particles and fluid are solved on an inertial time scale, which proves to be surprisingly simple and efficient, avoiding the costly linear algebra associated with Brownian dynamics. Complex confined geometries can be represented by a straightforward mapping of the boundary surfaces onto a regular three-dimensional grid. The hydrodynamic interactions between monomers are shown to compare well with solutions of the Stokes equations down to distances of the order of the grid spacing. Numerical results are presented for the radius of gyration, end-to-end distance, and diffusion coefficient of an isolated polymer chain, ranging from 16 to 1024 monomers in length. The simulations are in excellent agreement with renormalization group calculations for an excluded volume chain. We show that hydrodynamic interactions in large polymers can be systematically coarse-grained to substantially reduce the computational cost of the simulation. Finally, we examine the effects of confinement and flow on the polymer distribution and diffusion constant in a narrow channel. Our results support the qualitative conclusions of recent Brownian dynamics simulations of confined polymers [Jendrejack et al., J. Chem. Phys. 119, 1165 (2003) and Jendrejack et al., J. Chem. Phys. 120, 2513 (2004)].     

Simulations of polyelectrolyte dynamics in an externally applied electric field in confined geometry

We consider the dynamics of charged polymers in free solution in a slit geometry under the influence of an electrical field, applied at an angle to the plane parallel walls of the confinement. The simulations are carried out using the Brownian dynamics method with explicit counterions and implicit hydrodynamics. The hydrodynamic interactions between all the particles and the plane parallel walls are taken into account using a diffusion matrix which depends on slit geometry and the actual polyelectrolyte-solute conformations. We observe a selective transport of the charged polymers, as a function of the degree of polymerization and slit height.     

Anisotropic friction model of DNA electrophoresis in polymer solutions: Comparison with direct observations

The electrophoresis of DNA chains in uncrosslinked polymer solutions with a Brownian dynamics simulation with an anisotropic friction tensor was analyzed. According to the degree of anisotropy, three types of migration behavior are obtained: fluctuation without or with periodicity between U‐shaped and compact conformations, or migration with linear conformation. We found good agreement between our simulation results and the direct observations of DNA by fluorescence microscopy. © 2003 Wiley Periodicals, Inc. J Polym Sci Part B: Polym Phys 41: 1316–1322, 2003     

Dynamic Monte Carlo simulations of anisotropic colloids

We put forward a simple procedure for extracting dynamical information from Monte Carlo simulations, by appropriate matching of the short-time diffusion tensor with its infinite-dilution limit counterpart, which is supposed to be known. This approach - discarding hydrodynamics interactions - first allows us to improve the efficiency of previous dynamic Monte Carlo algorithms for spherical Brownian particles. In the second step, we address the case of anisotropic colloids with orientational degrees of freedom. As an illustration, we present a detailed study of the dynamics of thin platelets, with emphasis on long-time diffusion and orientational correlations.     

Calculation of diffusion coefficients for carbon dioxide + solute system near the critical conditions by non-equilibrium molecular dynamics simulation

A non-equilibrium molecular dynamics simulation was adopted to calculate the diffusion coefficients for a pseudo-binary system of carbon dioxide and for a carbon dioxide + solute system at 308.2 and 318.2 K. The calculated results were compared with the self- and tracer diffusion coefficients calculated by an equilibrium molecular dynamics simulation. The simulated results for the pseudo-binary system of carbon dioxide by the non-equilibrium molecular dynamics simulation are in good agreement with the results of self diffusion coefficients for pure carbon dioxide by the equilibrium molecular dynamics simulation. The simulated results of mutual diffusion coefficients for the carbon dioxide + solute system by the non-equilibrium molecular dynamics simulation are slightly lower than the results of the tracer diffusion coefficients by the equilibrium molecular dynamics simulation. The anomalous behavior of diffusion coefficients near the critical concentration was represented by the results of the non-equilibrium molecular dynamics simulation.     

Dynamics of polymers in a particle-based mesoscopic solvent

We study the dynamics of flexible polymer chains in solution by combining multiparticle-collision dynamics (MPCD), a mesoscale simulation method, and molecular-dynamics simulations. Polymers with and without excluded-volume interactions are considered. With an appropriate choice of the collision time step for the MPCD solvent, hydrodynamic interactions build up properly. For the center-of-mass diffusion coefficient, scaling with respect to polymer length is found to hold already for rather short chains. The center-of-mass velocity autocorrelation function displays a long-time tail which decays algebraically as (Dt)(-3/2) as a function of time t, where D is the diffusion coefficient. The analysis of the intramolecular dynamics in terms of Rouse modes yields excellent agreement between simulation data and results of the Zimm model for the mode-number dependence of the mode-amplitude correlation functions.     

Mass transport of O2 and N2 in nanoporous carbon (C168 schwarzite) using a quantum mechanical force field and molecular dynamics simulations

A hierarchical approach is used to calculate the single-component fluxes of N2 and O2 in nanoporous carbon molecular sieves (represented by C168 schwarzite) over a wide range of pressures and pressure drops. The self- and corrected diffusivities are calculated using equilibrium molecular dynamics simulations with force fields for the gas-carbon interactions obtained from quantum mechanical calculations. These results are combined with previously reported adsorption isotherms of N2 and O2 in C168 to obtain transport diffusivities and, by use of the Fick's equation of mass transport, to obtain single-component fluxes across the membrane. The diffusion coefficients and fluxes are also calculated using an empirical potential, which has been obtained by fitting low coverage adsorption data of N2 and O2 on a planar graphite sheet. By analyzing the diffusivities calculated with the ab initio potential in the limit of infinite dilution over the temperature range from 80 to 450 K, it is observed that the N2/O2 separation is energetically driven and a high selectivity of O2 over N2 can be obtained at low temperatures. However, with the empirical potential both the energetic and entropic contributions to selectivity were found to be close to unity. Similarly, by calculating single-component fluxes and ideal selectivities at 300 K and finite pressures it is found that the ab initio potential better explains the large O2/N2 selectivities of similarly sized molecules that have been observed experimentally. An interesting reversal in ideal selectivity is observed by adjusting the pressure at the two ends of the membrane. As a consequence, we predict that a highly selective kinetic separation in favor of either nitrogen or oxygen could be obtained with the same membrane depending on the operating conditions.     

Modeling gas permeation through membranes by kinetic Monte Carlo: applications to H2, O2, and N2 in hydrated Nafion®

We present a simulation tool in order to predict gas permeation through heterogeneous, microphase separated structures. The method combines dissipative particle dynamics (DPD) with kinetic Monte Carlo (KMC). Morphologies obtained from DPD are mapped onto a high density grid on which gas diffusion takes place. Required input parameters for the KMC calculations are the gas solubility and gas diffusion constant within each of the pure phase components. Our method was tested and validated for permeation of H(2), O(2), and N(2) gasses through hydrated Nafion membranes at various temperatures and water contents. We predict that membranes that contain an equal volume fraction of water, those with the highest ion exchange capacity exhibit the largest N(2) and O(2) permeation rates. For membranes of the same ion exchange capacity the H(2), O(2), and N(2) and permeability increases approximately linearly with Bragg spacing. We also predict that O(2) gas permeation depends much more on bottleneck phenomena within the phase separated morphologies than H(2) gas permeation. Overall, the calculated H(2) and O(2) permeability is found to be slightly lower than experimental values. This is attributed to the robustness of DPD resulting in ～7% larger Bragg spacing as compared with experiment and∕or increased gas solubility within the polymer phase with water uptake.     

Gas sorption and barrier properties of polymeric membranes from molecular dynamics and Monte Carlo simulations

It is important for many industrial processes to design new materials with improved selective permeability properties. Besides diffusion, the molecule's solubility contributes largely to the overall permeation process. This study presents a method to calculate solubility coefficients of gases such as O2, H2O (vapor), N2, and CO2 in polymeric matrices from simulation methods (Molecular Dynamics and Monte Carlo) using first principle predictions. The generation and equilibration (annealing) of five polymer models (polypropylene, polyvinyl alcohol, polyvinyl dichloride, polyvinyl chloride-trifluoroethylene, and polyethylene terephtalate) are extensively described. For each polymer, the average density and Hansen solubilities over a set of ten samples compare well with experimental data. For polyethylene terephtalate, the average properties between a small (n = 10) and a large (n = 100) set are compared. Boltzmann averages and probability density distributions of binding and strain energies indicate that the smaller set is biased in sampling configurations with higher energies. However, the sample with the lowest cohesive energy density from the smaller set is representative of the average of the larger set. Density-wise, low molecular weight polymers tend to have on average lower densities. Infinite molecular weight samples do however provide a very good representation of the experimental density. Solubility constants calculated with two ensembles (grand canonical and Henry's constant) are equivalent within 20%. For each polymer sample, the solubility constant is then calculated using the faster (10x) Henry's constant ensemble (HCE) from 150 ps of NPT dynamics of the polymer matrix. The influence of various factors (bad contact fraction, number of iterations) on the accuracy of Henry's constant is discussed. To validate the calculations against experimental results, the solubilities of nitrogen and carbon dioxide in polypropylene are examined over a range of temperatures between 250 and 650 K. The magnitudes of the calculated solubilities agree well with experimental results, and the trends with temperature are predicted correctly. The HCE method is used to predict the solubility constants at 298 K of water vapor and oxygen. The water vapor solubilities follow more closely the experimental trend of permeabilities, both ranging over 4 orders of magnitude. For oxygen, the calculated values do not follow entirely the experimental trend of permeabilities, most probably because at this temperature some of the polymers are in the glassy regime and thus are diffusion dominated. Our study also concludes large confidence limits are associated with the calculated Henry's constants. By investigating several factors (terminal ends of the polymer chains, void distribution, etc.), we conclude that the large confidence limits are intimately related to the polymer's conformational changes caused by thermal fluctuations and have to be regarded--at least at microscale--as a characteristic of each polymer and the nature of its interaction with the solute. Reducing the mobility of the polymer matrix as well as controlling the distribution of the free (occupiable) volume would act as mechanisms toward lowering both the gas solubility and the diffusion coefficients.     

Viscoelasticity of stretched polymer chains: Analytical theory and computer-aided simulation

The viscoelasticity of stretched polymer chains has been studied by the method of collisional dynamics. To this end, time correlation functions of the fluctuations of the microscopic stress tensor are modeled and relaxation moduli are expressed. Before, for stretched polymer networks, correlation functions used to be calculated in terms of an approximate theory that allowed one to estimate the strain dependences of loss modulus. The calculated dependences are shown to agree qualitatively with the results of measurements performed over a wide interval of strains, including prefracture strains. This theory is verified by comparing the time correlation functions of stress tensor fluctuations for a single stretched chain; these functions are found by computer-aided simulation and calculated on the basis of the existing analytical theory. In this case, a simple theory is adopted according to which a polymer molecule represents a chain composed of N atoms connected by freely jointed elastic bonds. The first and Nth atoms of this chain are attached by harmonic springs to immobile points located at a fixed distance. The decay of time correlation functions under study can be resolved into three stages. After a short initial interval provided by local motions, one can observe a region of power-law decay, which is followed by monoexponential decay at long times. The results of computer-aided simulation generally agree with the predictions of analytical theory. Certain discrepancies primarily concern the dependences of the exponent of power-law relaxation on the degree of chain stretching.