Conformationally averaged iterative Brownian dynamics simulations of semidilute polymer solutions Free

Polymer solution dynamics are of considerable relevance to both polymer processing applications¹ and fundamental polymer science.² Flows strongly affect the out-of-equilibrium conformational motions of polymer molecules and conversely these motions affect the rheological properties of the fluid. This interplay between the solvent and polymer has a significant impact on material properties in solution processing. Of particular interest is the semidilute unentangled concentration regime, which is characterized by polymer molecules that “overlap” and exhibit a large number of intermolecular interactions.^2,3 This contrasts from the dilute case, where individual polymers are isolated and primarily experience intra-molecular interactions. Processing applications of these semidilute solutions generally involve flows which take the solution significantly out of equilibrium,^4,5 but even in the equilibrium case, molecular diffusion and relaxation times depend strongly on the polymer concentration. Due to the ubiquity of semidiliute polymer solutions, it is important to have physical insights into the molecular origin of their dynamics.

To develop this insight, experimental techniques such as forced Rayleigh scattering,^6,7 pulsed field gradient nuclear magnetic resonance,^8,9 dynamic light scattering,^10,11 and fluorescence correlation spectroscopy¹² have all captured bulk solution dynamics. These measurements are commonly interpreted in the context of fundamental theoretical models that predict the molecular weight and concentration-dependence of the polymer diffusion coefficient.^2,3 Both scaling-based¹³ and kinetic theories^14,15 can capture the structure and dynamics of semidilute solutions; namely, they show how polymer solutions transition from Zimm behavior¹⁶ (used in dilute solutions) to Rouse¹⁷ and reptation-tube behaviors¹⁸ (used for concentrated solutions). This transition is conceptually captured via a length scale ξ, known as a correlation “blob” that is set by the polymer concentration. Zimm behavior occurs at length scales <ξ and concentrated (Rouse or reptation) behavior occurs for length scales >ξ. In equilibrium scaling arguments, this length scale is typically associated with both structural and hydrodynamic correlations, which are separate quantities but differ only by a factor of order unity. These ideas hold for asymptotic limits far from the crossover regions. Generalized scalings accounting for the double-crossover in solvent quality and concentration have also been explored by Brownian dynamics (BD) simulations.¹⁹ 

Bulk experimental measurements are only an indirect validation of these molecular theories. To directly connect polymer solution dynamics to molecular motions, fluorescently tagged DNA chains can be directly visualized in real time by single molecule fluorescence microscopy.²⁰ These single molecule experiments have provided evidence supporting theoretical models including reptation in melts,^21,22 stretching under transient extensional flows,^23–25 and tumbling in shear flows.^26–28 Most experiments have focused on the dilute and concentrated regimes, but recent single molecule experiments have investigated the transition from semidilute to concentrated as well.^29–32 These results are similarly in agreement with the existing scaling picture of semidilute solutions.^2,3 It remains unclear, however, how to extend these concepts to more complicated situations, such as semidilute polymers with non-linear architectures,³³ self-assembling structures,³⁴ or more generally semidilute solutions driven out-of-equilibrium. These are persistent challenges in understanding the coupling between polymer conformations and flow.

As an alternative to experiment and theory, simulations are a powerful tool to obtain molecular-level information of polymer dynamics. They allow for direct interrogation of the underlying theoretical picture with far fewer assumptions and resolve molecular degrees of freedom necessary to understand more complicated or out-of-equilibrium polymer solutions. Simulations are particularly valuable in the semidilute regime, where accurate modeling can provide further insight into polymer conformational dynamics. For example, simulations have found that the transition from Zimm to Rouse dynamics in non-dilute solutions is time-dependent as well as length-dependent.³⁵ Simulations of semidilute solutions in shear and extensional flows show that the effect of concentration can be accounted for by a rescaled Weissenburg number³⁶ which has since been observed in experiment.³² This effect has been confirmed by another simulation study that further explores the complex relationship between deformation, excluded volume, and hydrodynamics under shear flow.³⁷ Many predictions of scaling theory have been verified by simulation and connected to molecular forces as well.^35–39 

Simulations of semidilute polymer solutions remain limited, however, because even coarse-grained polymer models must include the effects of solvent dynamics for the long computation times needed to capture slow polymer relaxation modes. To capture both solvent and polymer time scales, an “effective” solvent is used in which a stochastic displacement accounts for solvent collisions. This solvent also must reproduce solvent-mediated interactions between particles, known as hydrodynamic interactions (HIs).² HI arises because particle motions introduce disturbance velocities in the solvent velocity field, affecting other particles.² The effective solvent is generally treated in one of two ways: (i) by quasi-explicitly resolving the solvent degrees of freedom and thus capturing HI by the exchange of momentum between the particles and solvent and (ii) by assuming that the solvent is a continuum and HI can be represented by an analytical expression for the disturbance velocity.² The advantage of the first approach is that the operation count of computing exchange between the particles and solvent is O(N).^40–44 The challenge in this class of simulations is to ensure the quasi-explicit solvent correctly resolves the fluid inertia by careful control of the time step. In particular, previous work has shown that accounting for vorticity diffusion at all length scales in overdamped systems leads to O(N^5/3) scaling.⁴⁵ 

The second approach, known as Brownian dynamics (BD), can be advantageous because it accounts for the fluid inertia implicitly. In this case, the Langevin equation motion for the position r_i of a polymer bead i at equilibrium is² 

where the first term accounts for pair-wise forces and HI via the diffusion tensor D. The second term accounts for the solvent collisions with the particle by a Gaussian random force f. Fluctuation dissipation requires that the noise is correlated to the HI as D = BB^T.² Therefore, BD simulations must involve the decomposition of the diffusion tensor. The decomposition is typically the bottleneck in BD simulations.⁴⁶ An exact implementation using Cholesky factorization leads to computational cost which scales with the number of particles as O(N³). More efficient techniques such as Fixman’s method^47–49 and the Krylov subspace⁴⁶ use polynomial approximations to take advantage of the fact that only the product of the decomposition and noise vector B·f is required and not the decomposition separately. These approximations reduce the scaling of the decomposition to O(N^2−2.25).^46,47

In semidilute BD simulations, another significant cost is the construction of the diffusion tensor itself. In this case, an Ewald sum which splits the slowly convergent HI into exponentially convergent real space and Fourier space contributions is necessary. An optimized pair-wise implementation leads to O(N^1.5) scaling of the Ewald sum or O(N^1.8) overall when Fixman’s method⁴⁷ is used to compute the Brownian noise. Clearly, these methods show superlinear scaling that limits BD simulations to several thousand particles. Alternatively, the HI can be computed by a particle mesh Ewald (PME) sum to achieve O(N log N) scaling.⁵⁰ This technique has been implemented in BD simulations which achieve similar scaling^39,51,52 via the Krylov subspace with PME. Ideas from explicit solvent simulations have also been implemented in BD to achieve approximately linear scaling with the number of particles.⁴⁵ Finally, alternative methods such as Lattice Boltzmann (LB) or multiparticle collision dynamics (MPCD) that use quasi-explicit solvent may be competitive, as shown by comparisons performed by Jain et al.³⁸ However, this comparison was performed prior to many recent developments,^39,45 so there is a need for further comparison. Nevertheless, this result suggests that there may still be significant opportunities to accelerate BD simulations.

We propose an alternative method of accelerating BD simulations based on physical rather than numerical approximations. Specifically, we take inspiration from the classic Zimm model,¹⁶ which pre-averages HI over the equilibrium chain conformations. Extensive effort has historically sought to extend this idea beyond isolated chains in equilibrium. For example, de Gennes proposed a model for pre-averaging the HI over the out-of-equilibrium conformation of the chain.⁵³ This approach led to the generalized Zimm model first considered by Öttinger⁵⁴ and later by Magda et al.⁵⁵ Öttinger further included Gaussian fluctuations in the model to account for perturbations of the chain conformation from the average,⁵⁶ and Prakash developed a more efficient approximation which allows for investigation of long chains and extrapolation to universal viscometric functions.^57,58 We have shown that this concept can be extended to Brownian dynamics simulations in an iterative Conformational Averaging (CA) method.⁵⁹ The method has been implemented for dilute linear polymer solutions at equilibrium and under extensional and shear flows and can be generalized to non-linear polymer architectures such as rings. The computational advantage of the CA method is that the Ewald sum and decomposition do not need to be performed repeatedly throughout the simulation. Instead, the average HI and decomposition are sampled over the course of the simulations and then averaged. An additional simulation is then performed using this average to obtain a self-consistent measure of the average HI in an iterative procedure outlined in detail in Sec. II. The sampling frequency is low enough that the simulation is computationally limited by the integration of the Langevin equation rather than the calculation of the HI and decomposition.

In this work, we extend the CA method to semidilute solutions, where fluctuations in the intermolecular HI preclude simply averaging the diffusion tensor over the polymer conformation as in the dilute case. Here, we introduce modifications that allow the intermolecular HI to vary over time while retaining a fast approximate measure of the HI. Additionally, we use the truncated expansion ansatz (TEA) to pre-average the decomposition and allow for fast calculation of the Brownian noise. The organization of this article is as follows: In Sec. II, we summarize the governing equations for BD simulations and introduce our iterative conformational averaging (CA) method. In Sec. III, we compare the computational cost and predictions of our method to those of other methods, analyze the results through static and dynamic properties, and compare to experimental results. We also investigate the screening of HI by the interaction of polymer molecules with a point force in solution.

We perform equilibrium BD simulations on polymer solutions of N_c chains and N_b coarse-grained beads per chain such that the total number of beads is N = N_bN_c. The beads have positions r_i, which move according to the Langevin equation

Tildes denote dimensionless quantities; positions are normalized by the bead radius (⁠$r̃=r/a$⁠), energies are normalized by the thermal energy k_BT [$Ũ$ = U/(k_BT)], times are normalized by the single-bead diffusion time [$t̃=t/τ0$⁠, where τ₀ = 6πηa³/(k_BT) and η is the solvent viscosity], and the mobility matrix is normalized by the drag coefficient of the spherical polymer beads [$D̃ij=Dij/(6πηa)$].

The beads interact via a potential $Ũ$ = $ŨS$ + $ŨB$ + $ŨEV$ that consists of connectivity, flexibility, and excluded volume, respectively. For connectivity, we consider a Hookean spring potential,

where $r̃ij=|r̃j−r̃i|$ is the distance between beads i and j and a large spring constant $κ̃S=200$ is chosen so that the beads are effectively connected by stiff springs. Flexibility is modeled by a bending potential

where θ_i is the angle formed between the bond vectors $r̃i,i−1$ and $r̃i+1,i$ and θ₀ = 0 is the equilibrium angle. The bending constant κ_B = 3.33 is chosen to give a persistence length of approximately 3 polymer beads. Excluded volume interactions are modeled by a Lennard-Jones (LJ) potential

The strength of the potential ũ = 0.1 is chosen so that the chain statistics are representative of a good solvent with the excluded volume scaling exponent ν ≈ 0.58. This value is determined from the results of single chain simulations using the scaling relation $⟨R̃g0⟩∼Nbν$⁠, where $⟨R̃g0⟩$ is the equilibrium radius of gyration in the dilute limit.

The solvent-mediated HI and Stokes drag are included by the mobility matrix given by the Rotne-Prager Yamakawa (RPY) tensor,^60,61

where $r^ij=r̃ij/rij$ is a unit vector in the direction of $r̃ij=r̃j−r̃i$ and I is the identity matrix. To satisfy the fluctuation-dissipation theorem, the Brownian noise $ξ̃i$ is related to the mobility tensor by $⟨ξ̃i(t)ξ̃j(t′)⟩=2D̃ijδ(t−t′)$ and $⟨ξ̃i(t)⟩=0$⁠. This requires the decomposition of the mobility matrix as $D̃=BBT$⁠. The Brownian noise is then $ξ̃i=2Bijfj$⁠, where f_j is a Gaussian random variable with mean 0 and variance 1.

Polymers are simulated in a cubic box of size $Ṽ=L̃3$ with periodic boundary conditions. The size of the box is determined by $Ṽ=N/c̃$⁠, where $c̃$ is the monomer concentration. We set the normalized concentration $c̃/c̃*$ where $c̃*=Nb/(4/3π⟨R̃g0⟩3)$ is the overlap concentration.

The long-range nature of HI is accounted for by an Ewald summed RPY tensor as formulated by Beenakker⁶² and later corrected for bead overlap by Jain et al.³⁸ This accounts for interactions between beads inside the primary image simulation box and all periodic replicas of the primary image. The Ewald sum overcomes the slow convergence of the HI by splitting the sum into real space and reciprocal space parts which both decay exponentially. Accounting for the correction due to self-interaction gives the mobility matrix in three terms

The parameter α determines the computational load of the real space and reciprocal space sums by their convergence rates. We choose $α=6/L̃$ for our simulations which has been shown to provide a good balance between the two sums.³⁸ The vector n = (n_x, n_y, n_z) with integer components specifies all images including the primary image. The prime on the sum over n in the real space term indicates that the primary image n = (0, 0, 0) is omitted for i = j. $Mα(1)$ is a 3 × 3 matrix and a function of the vector between the images of beads i and j, $r̃ij,n=r̃j−r̃i+nL̃$⁠,

The parameters C₁…C₄ are given by

In the reciprocal space sum, $k̃λ=2πL̃l$ and l ∈ Z³. The 3 × 3 matrix $Mα(2)$ is given as

where $k̃=|k̃|$ and $mα(2)$ is given as

The Ewald sum is optimized following the procedure of Jain et al.³⁸ to minimize computational time while ensuring the sum has converged to an error of approximately 10⁻⁴.

The computational bottleneck in BD simulations is typically the decomposition of the mobility matrix to calculate the correlated Brownian noise. Progress towards overcoming this bottleneck by mathematical approximations is discussed in the Introduction. In this work, we instead make a physical approximation by assuming that hydrodynamic interactions can be conformationally averaged (CA) over the structure of the polymer chain and the neighboring molecules. An averaged form of the Brownian noise is also used to maintain fluctuation-dissipation. The advantage of this approximation is that the HI and Brownian noise do not need to be explicitly computed throughout the simulation. The CA method requires relatively infrequent sampling of the HI to provide a sufficient average over the conformational space of the solution structure. The challenge is that there is a self-consistent relationship between the average HI and the solution structure such that the average HI generally cannot be determined before simulation. We overcome this by performing iterative BD simulations until converged to a self-consistent result. Previously, we have implemented this iterative procedure for single chain BD simulations using an approximation of the Langevin equation

This differs from the original Langevin equation in that we have introduced a mobility matrix and a Brownian noise matrix which are averaged over polymer conformations. Fluctuation-dissipation is maintained by the new random velocity $ξ̃i(w)$⁠, which satisfies $⟨ξ̃i(w)(t)ξ̃j(w)(t′)⟩=2⟨D̃ij⟩(w)δ(t−t′)$⁠. In addition, a superscript (⁠$w$⁠) has been introduced to denote the iteration number, with the first simulation as $w$ = 0. For a complete description of the single chain method, we reference our previous publication;⁵⁹ here, we briefly outline the method:

1. Start with a guess for the averaged mobility matrix, typically the freely draining (FD) case, $⟨D̃ij⟩(0)=δijI$⁠. It is also possible to make guesses which are closer to the converged average, such as the predictions of Zimm theory.

2. Perform a BD simulation using the averaged mobility matrix to capture the solvent-mediated forces and Brownian noise. The decomposition of the mobility matrix is only required once at the start of each iteration.

3. Use the BD simulation to calculate the time averaged mobility matrix for this iteration. Assuming steady state, the phase space average can be replaced with a time average where the conformational probability is discretely sampled via the BD simulation.

4. Repeat steps 2 and 3 with increasing values of $w$ = 1, 2, 3, … until converged to a self-consistent mobility matrix $⟨D̃ij⟩(w+1)≈⟨D̃ij⟩(w)$⁠.

Each iteration is run for 15τ_R, where τ_R is the longest relaxation time of a single chain. We find this provides sufficient sampling of conformational space. In the first iteration using freely draining dynamics, τ_R is the Rouse relaxation time. For subsequent iterations including HI, this is the Zimm relaxation time.

Prior work by the authors has shown that the CA method quantitatively matches standard BD simulations for single linear and ring polymers in equilibrium and out-of-equilibrium under steady state flow conditions with the exception of strong shear flows that exhibit large conformational oscillations.⁵⁹ Semidilute solutions are more challenging, however, because CA based on the 3N × 3N bead-index based diffusion tensor would remove any relative spatial information for a pair of two hydrodynamically interacting monomers on different chains. Instead, we introduce a spatial average to capture intermolecular HI, while retaining the previous index-based procedure for intramolecular HI. To modify our approach, we introduce notation to specify both the bead indices i and j (running from i, j = 1 to i, j = N_b) and chain indices α and β (running from α, β = 1 to α, β = N_c). The mobility matrix is thus fully denoted by the bead and chain indices of interacting particles, their spatial separation, and the iteration number $w$⁠; to capture all of these attributes, we write the full mobility as $D̃ijαβ,(w)(r)$⁠, which is given by the Rotne-Prager-Yamaka tensor discussed previously. The self-HI for individual chains is calculated in the same fashion as in prior work by the authors,

Here, the superscript S denotes that this is a self-HI term, and while the average is over all chains, the Kronecker delta ensures that α = β. The denominator calculates the number of chains. The integral in the numerator is over all possible conformations of all chains, as denoted by the path integral. The value P^(⁠$w$⁠)({r^α}) is a configurational probability as described previously. In practice, we replace this conformational path integral with an integral over a steady-state trajectory,

This self-HI mobility $⟨D̃ijS,(w)⟩$ retains the indices i and j. By contrast, the intermolecular HI will retain the spatial location by averaging over all indices at a specific location r,

The superscript O indicates the intermolecular HI. The indices i, j, α, and β are all part of the sums, while a Dirac delta function considers only the bead pairs that contribute to the diffusivity at a separation of r. In the infinitesimal limit, $⟨D̃O(r)⟩(w)$ is simply the RPY tensor. However, we assign similar values of r to a discrete grid (Fig. 1). This transforms the mobility matrix from a function of the continuous distance between beads, $⟨D̃O(r̃)⟩(w)$⁠, to a function of the discrete grid space displacement, $⟨D̃O(Δr̃)⟩(w)$⁠, where $Δr̃=(Δx̃,Δỹ,Δz̃)$ is the displacement between beads rounded to the nearest grid point. In the subsequent iterations $w$ + 1, HI between chains is accounted for using the grid average measured in the previous iteration $w$⁠. At any given time, the displacements of the N(N − N_b)/2 bead pairs are calculated and rounded to points on the discrete grid. The resulting mobility matrix retains spatial information because the choice of the matrix element to apply depends on relative position $Δr̃ij$⁠; because we update which matrix element is used over the course of the simulation, this method incorporates fluctuations into the intermolecular portion of the HI. The resulting form of the Langevin equation gives the total velocity as a combination of the intramolecular and intermolecular components. The time-dependence of the intermolecular HI has been shown explicitly to reflect the change in $Δr̃ij$ with the conformations of surrounding chains,

Here we have constructed the 3N × 3N mobility matrix $⟨D̃eff⟩(w)$ from the two averages

The matrix multiplication of the mobility tensor with the bead forces $⟨D̃ijαβ,eff⟩∇r̃jU$ is then the same as in traditional BD simulations.

The quantitative details of the grid are chosen to minimize computational time while retaining an accurate measure of the continuous HI. We use a cubic grid with points separated by one bead diameter d_g = 2.0. The number of grid points in each dimension is the size of the box divided by the bead diameter plus an additional grid point allocated for the origin, $Md=L̃/dg+1$⁠. The total number of grid points is $M=Md3$ for the cubic simulation box. Out-of-equilibrium flow conditions may require modifications to the grid to account for deformed box boundary conditions, but in this work we only consider equilibrium dynamics. Large simulation boxes will have memory limitations because the number of grid points scales as $M∼L̃3$⁠. To address this, we use a fine grid at small displacements and a coarser grid at large displacements, d_g = 4.0, in simulations with large box sizes. A coarse grid is reasonable because for large displacements, the HI is slowly varying. Furthermore, we show in Sec. III E that the HI is often screened at these distances. The HI is sampled approximately 100 times per τ_R. This sampling interval is chosen so that the computationally expensive Ewald sum is not the bottleneck in the simulation while retaining an accurate average of conformational space. In the second and subsequent iterations, the fluctuating grid space mobility matrix is updated once every 10 time steps with a time step of 10⁻³τ₀. Similarly, this interval is long enough that the matrix multiplication of the mobility matrix with the pair-wise forces $⟨D̃ijαβ,eff⟩∇r̃jU$ remains the bottleneck while retaining an accurate measure of the HI. Previous work has shown that this interval could be increased to 50 or 100 time steps without significantly affecting dynamic properties.⁴⁶ 

We also need to account for fluctuations in the average Brownian noise. In single chain CA, the mobility matrix and Brownian noise matrix do not fluctuate, so we only need to carry out the Cholesky decomposition once per iteration. In the semidilute case, however, we retain fluctuations in the intermolecular HI; this would be a bottleneck similar to traditional BD simulations. To overcome this, we use the truncated expansion ansatz (TEA) to measure a time-averaged decomposition. The “ansatz” of the TEA is that the square root of the mobility matrix can be factorized into individual two-body contributions such that the correlated random vector in the Brownian noise takes the form

Following the notation of Geyer and Winter,⁶³ Eq. (19) gives one component of the size 3N Brownian noise vector. We use indices l and m to emphasize that this sum is over entries in a single row l of the mobility matrix and that Eq. (19) does not follow index notation. This contrasts with the bead index notation used previously, which denotes a 3 × 3 tensor relating the vector quantities of force and velocity for simulation beads. There are thus 3N coefficients C_l which ensure each bead recovers the correct self-mobility. The weighting factor β′ is chosen so that the Brownian noise approximates the mobility matrix decomposition by satisfying fluctuation-dissipation via $cov(ξ̃,ξ̃)≈D̃$⁠. For diagonal entries of the sum, corresponding to l = m, the weighting factor β′ = 1. All off diagonal entries have the same weighting factor, given by

where ϵ is an average over the off-diagonal entries of the mobility matrix

The coefficients are given by

Details on the derivation of the TEA can be found in the original paper by Geyer and Winter.⁶³ The conformationally averaged TEA is implemented by calculating the parameters β′ and C_l each time the Ewald sum mobility matrix is sampled. We thus calculate the time-average TEA parameters for a given iteration

Here, time-averages are over T samples of the mobility matrix during a steady-state simulation. We motivate the use of these averages instead of their instantaneous values by observing that that these values do not deviate significantly from these averages even when they were calculated at every time step. These quantities can then be used to calculate the Brownian noise as

The advantage of using the TEA is that once the averaged parameters are obtained, the calculation of the Brownian noise requires only the sum in Eq. (25). Using other methods, the expense would be similar to the full BD simulation because the decomposition would have to be repeated in full each time the mobility matrix is updated. The sum in the TEA can be represented as a matrix-vector product by first dividing the diagonal components of the mobility matrix by β′ such that β_ii = 1 to recover Stokes drag. Then the simulation implementation is

Furthermore, we can implement a “block” version of the TEA similar to the block TEA by using a block of s random vectors for f_m. Then the matrix-vector product becomes a matrix-matrix product between the 3N × 3N mobility matrix and 3N × s block random vectors to yield s Brownian noise vectors. This is advantageous because the mobility matrix generally changes slowly such that it does not need to be updated every time step. Then the update interval λ_RPY > 1 and we can use s = λ_RPY. The matrix-matrix product is faster than performing multiple matrix-vector products and has efficient implementations which can be significantly accelerated by multi-threading.

Simulations using the TEA are known to reproduce the dynamic properties of the Cholesky decomposition for dilute theta polymers.^46,63,64 The main assumption of the TEA is that HI coupling is weak and is therefore unsuitable for dense systems such as collapsed polymers and dense colloidal suspensions. We maintain a low volume fraction by simulating semiflexible polymers which have a larger radius of gyration and thus a lower overlap concentration. This allows us to reach high normalized concentrations c/c^* while remaining below a volume fraction of ϕ ≈ 0.10, above which the TEA begins to deviate significantly from the Cholesky decomposition.⁴⁶ We have tested the accuracy of the TEA with an error criterion⁴⁶ 

where $C=cov(ξ̃,ξ̃)$ is the covariance matrix of the Brownian noise and |·| is the Frobenius norm, defined as the square root of the sum of the squares of matrix entries. To follow fluctuation-dissipation, the covariance must satisfy $C≈D̃$⁠. We find that for the volume fractions used in this work, the error is close to that of the Krylov subspace method. Using the conformationally averaged TEA parameters gives an error criterion close to the original TEA. For higher density systems such as flexible, short polymers in semidilute theta solutions, the weak HI coupling assumption of the TEA may not be appropriate. In this case, the decomposition could be calculated instantaneously from the average HI using Fixman’s approximation or the Krylov subspace at an increased computational cost. In this work, we consider only the TEA.

We now compare the computational time of the CA method to the TEA and the Krylov subspace as a function of the number of polymer beads for N = 1000 − 10 000. For the CA method, the cost of the FD iteration and HI iteration is shown separately. The FD iteration scales approximately linearly with N and is an order of magnitude or more faster than the HI iteration. Therefore, the expense of conformationally averaging the HI in the initial FD simulation is relatively small and can be neglected when comparing the CA method to other techniques for moderate system sizes. The subsequent HI iteration cost scales as O(N²) and is an order of magnitude faster than the TEA or Krylov subspace simulations. The O(N²) scaling is due to the matrix-vector multiplication of the mobility matrix with the pair-wise forces. The reduction in computational time compared to the TEA and Krylov subspace simulations is due to the reduced frequency of performing the Ewald sum (typically every λ_RPY = 25 − 100 time steps in the TEA and Krylov methods versus every ∼10 000 time steps in the CA method) and the decomposition to find the TEA parameters β and C_i, which is performed at the same frequency as the Ewald sum.

Figure 2 shows the serially executed (1 thread) scaling with N for chains of N_b = 100 at a concentration c/c^* = 1.0 with λ_RPY = 25 and number of chains increasing from N_c = 10 − 100. The Krylov subspace simulation uses the block Krylov with error tolerance E_k < 0.01, which has been shown to be sufficient for BD simulations of polymer solutions.⁴⁶ The inset shows scaling of the computational time with the number of threads used. The Ewald sum has been parallelized using OpenMP. The matrix-vector product of the mobility matrix with the random vector and the matrix-matrix products in the TEA and block Krylov calculation of the Brownian noise use Intel Math Kernel Library (MKL) routines. The simulations are performed on dual 8-core Intel Xeon E5-2667 processors.

For larger systems, the O(N²) scaling causes the computational time for all methods compared here to become prohibitive. In this case, a particle-mesh Ewald technique is more appropriate due to the O(N log N) scaling. Comparing to the matrix-free technique of Liu and Chow⁵² implemented for polymer solutions by Sadaat and Khomami,³⁹ we observe a crossover at N ≈ 10⁵. The matrix-free technique considers a different model for the polymer interactions, so this comparison is not quantitative and is not shown in Fig. 2. However, it does indicate that the CA method significantly reduces the prefactor in the computational cost. Jain et al.³⁸ have shown that explicit solvent models such as LB and Stochastic Rotation Dynamics (SRD) are significantly faster than BD simulations and have O(N) scaling. BD simulations implementing PME have shown that the issue of scaling can largely be overcome.⁴⁵ Implementation of a PME technique in the CA method could reduce the O(N²) scaling while retaining the low prefactor, which would make BD simulations of polymer solutions comparable to explicit solvent models. Further optimizations for systems with O(10⁶) polymer beads are an interesting area of further study.

To validate the approximations of the CA method, we compare to non-CA BD simulations where the mobility matrix and its decomposition are recalculated every 25 time steps. We use the TEA and the Krylov subspace to rapidly perform this calculation and use a combination of static and dynamic properties to evaluate the accuracy of the CA procedure. For the former, we focus on the radius of gyration, defined for each polymer chain as $⟨Rg02⟩=1/N∑iNbri−rc2∼Nb2ν$ where r_c is the center of mass of the chain. The first equation provides a molecular definition of $⟨Rg02⟩$⁠, while the second relationship expresses the scaling prediction as a function of N_b and the Flory exponent ν. We capture the dynamic behavior of the polymer chains by calculating the long time diffusion constant defined as

and the zero-shear viscosity calculated via the Green-Kubo formula² 

where the stress tensor is determined by the Kirkwood forumula² 

We simulate a system with N_b = 50 and N_c = 20, which is small enough to be computationally tractable by both the TEA method and the Krylov subspace. We vary the concentration from 0.01 c/c^* to 3.0 c/c^* in order to remain at a low volume fraction for short chains with a high overlap concentration. More details on the simulation parameters can be found in Table I.

Figure 3 plots the radius of gyration normalized by the dilute limit value as a function of concentration. The radius of gyration in the FD and HI iterations of the CA method is nearly constant for all concentrations. Both the HI iteration and the TEA simulation predict the size of the chain to be slightly larger than in the FD simulation. The three simulations should agree if fluctuation-dissipation is strictly maintained, but there is a constant quantitative difference at all concentrations due to the TEA approximation that the Brownian noise can be constructed by pairwise contributions. The deviation in the radius of gyration does not qualitatively change the dynamics of the system, which have been shown to match well with standard BD simulations.^46,63,64 However, we again note that this method becomes less accurate for more concentrated simulation systems. The Krylov subspace simulations are in close agreement with the FD simulation. Additionally, the size of the chain is slightly larger in the CA simulations with HI than either the TEA or the Krylov simulations. We attribute this to the conformationally averaged intramolecular HI which uses a constant diffusion tensor for HI between polymer beads on the same chain. At high concentrations, the simulations are more sensitive to approximations such that neglecting intramolecular fluctuations causes a small change in the size of the chain.

We can also compare dynamic property predictions of the TEA, Krylov, and CA simulations. Figure 4 plots the center-of-mass, long-time diffusion constant D_L as a function of c/c^*. For the first iteration $w$ = 0, the freely draining mobility matrix results in slow diffusion (red circles) that weakly decreases above c/c^* ∼ 1. The subsequent iteration $w$ = 1 includes HI and there is a concomitant increase in D_L due to the presence hydrodynamic shielding. Each bead no longer experiences Stokes drag independently, as in the FD case. In this case, the value of D_L from both CA and TEA is in agreement (Fig. 4, black squares and blue triangles). Increasing concentration leads to a marked decrease in D_L, especially as c/c^* approaches the overlap concentration and enters the semidilute regime. An additional HI iteration, $w$ = 2, reveals that the dynamic properties have nearly converged after one iteration. This is expected; at equilibrium, chain conformations do not depend on hydrodynamics, so averaging the HI over an FD simulation is sufficient. There is a small quantitative change originating from the structural differences caused by the TEA as observed in the static properties.

The Krylov simulations again show quantitative differences from the TEA and CA simulations. At low concentrations, the diffusion constant from the Krylov simulations is approximately 5% larger, in agreement with previous work.^46,63,64 In the semidilute regime, the Krylov and TEA simulations agree and the CA method has a slightly higher diffusion constant. We also attribute this to the conformational averaging of the intramolecular HI as in the case of the radius of gyration. Fluctuations in the TEA parameters β′ and C_i are small, so deviations are caused primarily by the diffusion tensor itself and not the decomposition technique.

To further investigate the effect of fluctuations on the dynamics, we also consider the zero-shear viscosity. The diffusion constant is relatively insensitive to preaveraging as shown by the success of Zimm theory, but Öttinger has found that incorporating Gaussian fluctuations into Zimm theory leads to a decrease in the viscosity.⁵⁶ In Fig. 5, we find a similar result. All simulations show similar scaling with concentration, but the freely draining iteration of the CA method shows the highest viscosity as expected because the polymer is unshielded from solvent drag. Simulations including HI have a lower viscosity, with the Krylov simulations being slightly less than the CA method. While we do not draw any quantitative conclusions due to the difficulty in measuring the long-time terminal regime of the stress-stress autocorrelation function used in the Green-Kubo relation, we speculate that this difference reflects the absence of fluctuations in the intra-molecular portion of the diffusion tensor.

We now consider further predictions of the method in comparison to experiment and theory. We simulate chains of length N_b = 100 and N_b = 150 and compare to the results for N_b = 50 from Sec. III B to study the effect of molecular weight. Detailed simulation parameters and inputs can be found in Table I. For all but the most concentrated simulations, the box size is chosen such that the average polymer end-to-end length is less than the minimum image distance $⟨R⋅R⟩<L̃/2$⁠. We also emphasize that the solutions considered here are unentangled because the chains are too short and/or dilute. The concentration above which polymers begin to entangle c_e is typically in the range of 4–30 c/c^*,⁶⁵ but the highest concentration simulated here is 5.77 c/c^*. We again consider the radius of gyration as a function of concentration. In the standard scaling picture for flexible chains, the chain conformation is determined by solvent statistics associated with excluded volume below a correlation length ξ_C. Above the correlation length, excluded volume interactions are screened and the polymer again follows random walk behavior leading to the prediction⁶⁵ $⟨Rg⟩=⟨Rg,0⟩c/c*−(2ν−1)/(6ν−2)$⁠. For a good solvent with ν ≈ 0.588, this gives $⟨Rg⟩∼(c/c*)−0.12$⁠. In contrast to this prediction, which has been observed in prior simulations of flexible polymers,³⁷ we are simulating semiflexible polymers with a relatively high persistence length, so this scaling is not observed. Instead, we seem to be in the so-called “marginal solvent regime” where the persistence length is sufficiently long that scaling theories break down and dynamic properties show mean-field behavior.⁶⁶ Mean-field theory predicts no concentration dependence, or $⟨Rg2⟩∼(c/c*)0$⁠.⁶⁶ This is close to the behavior we observe in simulation (Fig. 6). For longer chains, we begin to see a weak concentration dependence associated with the transition from good to marginal solvents for the larger thermal blob size. Again, we see a small quantitative difference between the FD and HI iterations associated with the TEA that does not affect the concentration dependence.

Scaling theory can also be used to predict dynamic properties. The analysis is similar, with the chain following dilute behavior below the correlation length ξ_C and concentrated behavior above ξ_C. This is based on the expectation that the correlation length and the hydrodynamic correlation length are proportional, ξ_C ∼ ξ_h. In this case, unscreened hydrodynamics (i.e., Zimm dynamics) are observed for length scales <ξ_h, while HI is screened for length scales >ξ_h (i.e., Rouse dynamics). These results lead to the scaling relationship⁶⁵ 

Here, $DL∼(c/c*)−0.54$ for good solvents. This contrasts with the behavior of D_L expected in the absence of HI, which should be independent of concentration (⁠$DL∼(c/c*)0$⁠) because the Stokes drag on each monomer does not depend on its environment, except when monomer-monomer excluded volume constraints become important. In Fig. 7, we plot the diffusion constant normalized by the dilute limit value from single chain simulations (D_L/D₀). We observe that for $w$ = 0 (FD), there is only a weak dependence of D_L on c/c^*, consistent with the expected $(c/c*)0$⁠, scaling. At large values of c/c^*, however, we find that there is a measurable decrease in D_L for all N_b. We attribute this to excluded volume, which should be the largest for N_b = 50 where the volume fraction of beads is the largest; indeed, this case exhibits the largest decrease in D_L.

When HI is included (⁠$w$ = 1), there is a pronounced decrease in the value of D_L, reflecting the change in the HI environment. Similar to the FD case, we observe an enhanced decrease in D_L compared to the longer polymers for the short, N_b = 50 chains. This is again consistent with our hypothesis that excluded volume plays a large role in dynamics. In agreement with scaling predictions, longer chains (N_b = 100, 150) collapse onto a single curve when D_L/D₀ is plotted versus (c/c^*). We note that we obtain the scaling $DL∼(c/c*)−0.6$⁠, which Eq. (31) predicts corresponds to a value of ν = 0.57. This is close to the prediction for a good solvent; however, it is unclear how this is modified by (1) the marginal solvent behavior observed in equilibrium and (2) the excluded volume effects seen for FD (⁠$w$ = 0).

We find that inclusion of HI allows us to capture the diffusion behavior observed in experiment. Measurements from single molecule DNA experiments³⁰ and dynamic light scattering of polystyrene¹¹ are co-plotted with the simulation results. The FD simulation does not show the correct concentration dependence. The simulation with conformationally averaged HI shows the same scaling as the experimental data and quantitatively matches the single molecule experiments (blue upside-down triangles). The light scattering data shows the same scaling but an earlier transition into the semidilute regime.¹¹ This could be related to differences (1) due to the flexibility of individual polymer chains used in the scattering experiment versus the semiflexible chains used in simulation and (2) in how the overlap concentration c^* is defined for light scattering experiments versus for our simulations.

Above we introduced HI screening in the context of scaling theory. Concentration screening is a fundamental concept in polymer physics and is widely utilized to describe the conformation of a polymer in semidilute solutions and melts. HI screening is less well studied but equally important for dynamic properties below c^**, the concentration where the correlation length shrinks below the thermal blob size. Calculations using scaling theory¹³ or effective medium theory^14,15 predict that the concentration screening length ξ_C and the HI screening length ξ_h have the same scaling with concentration. This is intuitive because they are both driven by the increase in monomer density. The scaling theory prediction for the concentration correlation length is

or $ξC∼(c/c*)−0.76$ for good solvents and $ξC∼(c/c*)−1.0$ for theta solvents.⁶⁵ Effective medium theory is consistent with this scaling theory for both theta¹⁵ and good¹⁴ solvents. We can directly and quantitatively measure HI screening and show that it behaves similar to theoretical predictions. This has important implications for the development of constitutive models for semidilute polymer solutions, where the effects of HI are generally accounted for by an effective friction factor.⁶⁷ 

In simulation, we measure HI screening by placing a point force in solution. The point force exerts a disturbance velocity on the surrounding polymers and solvent. We can then calculate the total disturbance velocity due to the point force and polymers at a displacement r by

where F^pf refers to the strength of the point force and $D̃pf$ refers to the 3 × 3 diffusion tensor at a displacement r. The first term gives the disturbance velocity in pure solvent, and the second term gives the effect of the background polymer solution with the diffusion tensor from CA, $⟨D̃O(Δr̃j*)⟩(w>0)$⁠. Here, the diffusion tensor is evaluated at $Δr̃j*$⁠, which is the discretized distance between the polymer bead j at a position $r̃j$ and a displacement $r̃$ from the point force. With this measurement, we can determine the HI screening length from the effective medium theory result¹⁴ 

where $v$_s is the pure solvent disturbance velocity and r = |r|. We can then use the pure solvent velocity from the RPY tensor and fit the theoretical prediction to the simulation result via the screening length ξ_h. We normalize the simulation results so that the disturbance velocity $v$_s(r = 0) at the point force location is given by the product of Stokes drag and the point force strength in agreement with the theoretical framework. We choose the magnitude of the point force so that it does not significantly perturb chain conformations out-of-equilibrium, yet is large enough that the resulting velocity perturbation is easily measurable.

The disturbance velocity also depends on the direction relative to the direction of the applied force as reflected in the RPY tensor. We choose to apply the force in the x-direction with the point force located at the origin. We then measure the disturbance velocity along the x-axis and y-axis. As seen in Figs. 8(a) and 8(b), the decay in the y-direction is faster than in the x-direction which is expected due to the form of the RPY tensor. This anisotropy is not reflected in most theoretical models, however, because they average over all directions of the perturbing force F^pf to reflect the isotropy of the polymer solution. In the equilibrium case, this feature does not appear to be essential as the scaling theory results agree with experiment. This will change when flows drive the system out-of-equilibrium, however, and the underlying solution is no longer anisotropic. This coupling between fluid anisotropy and the resulting HI is an intriguing topic for future study.

Fitting the screening length to the simulation results, we find that the screening length in both the x and y directions scales with concentration as $ξh∼(c/c*)−1$ [Fig. 8(c)]. This is in agreement with the effective medium theory result for the scaling of ⟨R_g⟩, which predicts a mean-field structure of the semidilute solution despite the good solvent due to the semiflexibility of the polymers. At concentrations lower than the overlap, the simulation results diverge as HI becomes unscreened in the dilute solution and the theoretical description of the effective medium becomes less accurate. At the overlap concentration, the correlation length is approximately equal to the radius of gyration, $ξC=⟨Rg02⟩1/2$⁠. This differs somewhat from the simulation results, although this could be due to the method used here to obtain ξ_h. Furthermore a small quantitative difference between the two length scales is expected.¹³ 

In this work, an iterative conformational averaging method for computing hydrodynamic interactions and the Brownian noise in polymer solutions has been proposed based upon the polymer configuration and the locations of surrounding molecules. This generalization of the dilute CA method⁵⁹ retains many of the same advantages. Specifically, the execution time per time step is reduced by an order of magnitude or more because the computationally expensive Ewald sum and decomposition are replaced by memory calls to the spatially averaged HI and TEA parameters, respectively. We have shown that this approach quantitatively matches the predictions of a traditional BD simulation and experiment. Additionally, we have used the method to investigate the concentration-dependence of HI screening in polymer solutions.

There are also limitations to the method that must be considered for the specific use case. While the dynamic properties generally match those of traditional BD simulations, we have found that neglecting fluctuations in the intramolecular HI can cause a quantitative deviation in the zero-shear viscosity. Additionally, static properties will deviate slightly due to the assumptions of the TEA. If highly accurate simulations are required, the resolution of the grid-space average HI is also important. In particular, the preliminary results from semidilute CA simulations in extensional flow suggest that polymer stretch and solution stress are sensitive to intramolecular fluctuations and the resolution of the grid-space HI. A complete analysis of these approximations and their effect on out-of-equilibrium dynamics will be carried out in a future publication.

This method is particularly relevant for the determination of HI screening associated with microstructural changes. Here we have considered only concentration and degree of polymerization, but topological effects relevant in branched and ring polymers are also of significance to industrial and biological research. For example, dilute ring polymer solutions show a delayed coil-stretch transition due enhanced hydrodynamic backflow associated with the looped conformation.⁶⁸ The concentration dependence of this phenomenon due to topological interactions with surrounding rings is an interesting area for further study. More generally, the idea of flow-dependent HI screening and generalization of the correlation length to account for flow-driven anisotropies of the polymer conformation is of importance for developing accurate constitutive models of polymer solutions.

Areas of further improvement for the CA method include reducing the O(N²) scaling and implementing a more accurate decomposition method. Currently, large scale simulation remains challenging relative to methods which achieve reduced scaling by FFT approximations to the wave space mobility. It is possible that a similar implementation or other algorithms such as the fast multipole method for the RPY tensor⁶⁹ could be extended to the CA method. The challenge would be to maintain the small prefactor that makes the CA method advantageous while achieving lower complexity. Integrating ideas from these methods into the CA algorithm could also improve the accuracy of the grid space intermolecular average. The use of the TEA for approximating the decomposition is also a concern for dense systems. Semidilute polymer solutions are characterized by low segmental density making the TEA an appropriate choice. For specific problems such as entangled solutions of relatively short chains, however, the assumption of weak HI may lead to more significant errors in the correlation between HI and Brownian noise. In this case, an average approximation to the eigenspectrum based upon iterative techniques may be required.

This work was funded by the National Science Foundation under Grant No. CBET-1803757. The authors also acknowledge helpful comments on the manuscript from Charles M. Schroeder.