Bead-bead interaction parameters in dissipative particle dynamics: relation to bead-size, solubility parameter, and surface tension

Dissipative particle dynamics (DPD) is a mesoscale modeling method for simulating equilibrium and dynamical properties of polymers in solution. The basic idea has been around for several decades in the form of bead-spring models. A few years ago, Groot and Warren established an important link between DPD and the Flory-Huggins chi-parameter theory for polymer solutions. We revisit the Groot-Warren theory and investigate the DPD interaction parameters as a function of bead size. In particular, we show a consistent scheme of computing the interfacial tension in a segregated binary mixture. Results for three systems chosen for illustration are in excellent agreement with experimental results. This opens the door for determining DPD interactions using interfacial tension as a fitting parameter.     

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)].     

On Coarse-Graining by the Inverse Monte Carlo Method: Dissipative Particle Dynamics Simulations Made to a Precise Tool in Soft Matter Modeling

We present a promising coarse-graining strategy for linking micro- and mesoscales of soft matter systems. The approach is based on effective pairwise interaction potentials obtained from detailed atomistic molecular dynamics (MD) simulations, which are then used in coarse-grained dissipative particle dynamics (DPD) simulations. Here, the effective potentials were obtained by applying the inverse Monte Carlo method [Lyubartsev and Laaksonen, Phys. Rev. E. 52, 3730 (1995)] on a chosen subset of degrees of freedom described in terms of radial distribution functions. In our first application of the method, the effective potentials were used in DPD simulations of aqueous NaCl solutions. With the same computational effort we were able to simulate systems of one order of magnitude larger than the MD simulations. The results from the MD and DPD simulations are in excellent agreement.     

Noncontinuum effects on the mobility of nanoparticles in unentangled polymer solutions

The results for the diffusivity of nanoparticles in unentangled semidilute polymer solutions obtained using coarse‐grained simulations are presented. The results indicate that for particle sizes smaller than the polymer radius of gyration, the nanoparticle diffusivities deviate from Stokes–Einstein predictions and depend explicitly on the polymer radius of gyration and the polymer solution correlation lengths. Scaling ideas proposed are invoked for rationalizing such noncontinuum effects and demonstrate that the simulation results could be collapsed onto a single universal function of the depletion thickness, the polymer radius of gyration, and the particle radius. © 2016 Wiley Periodicals, Inc. J. Polym. Sci., Part B: Polym. Phys. 2016 , 54, 2145–2150.     

Multiple "time step" Monte Carlo simulations: application to charged systems with Ewald summation

Recently, we have proposed an efficient scheme for Monte Carlo simulations, the multiple "time step" Monte Carlo (MTS-MC) [J. Chem. Phys. 117, 8203 (2002)] based on the separation of the potential interactions into two additive parts. In this paper, the structural and thermodynamic properties of the simple point charge water model combined with the Ewald sum are compared for the MTS-MC real-/reciprocal-space split of the Ewald summation and the common Metropolis Monte Carlo method. We report a number of observables as a function of CPU time calculated using MC and MTS-MC. The correlation functions indicate that speedups on the order of 4.5-7.5 can be obtained for systems of 108-500 waters for n=10 splitting parameter.     

Binding of poly(amido amine) dendrimer to sodium hyaluronate in aqueous NaCl solution

Binding of poly(amido amine) dendrimer to sodium hyaluronate (NaHA) in aqueous 0.25 M NaCl solution has been investigated by static light scattering. It was observed that the apparent weight-average molecular weight and the radius of gyration increase with the ratio of NH(2) terminal groups in the dendrimer to the carboxylate groups in the NaHA, [NH(2)]/[COO(-)]. Up to [NH(2)]/[COO(-)]=31, the observed variation of molecular weight was reproduced by the "average binding" model, where an average number of dendrimers binds to each NaHA chain. Based on the "critical binding" model, the maximum number, n(max), of dendrimers which can bind to a NaHA chain was calculated to be n(max)=300 for a solution of [NH(2)]/[COO(-)]=56. The obtained value corresponds to the binding of one dendrimer per 1.5 repeating units on a NaHA chain. It is suggested from the observed radius of gyration that, while the dendrimer-NaHA complexes of [NH(2)]/[COO(-)] up to 5 maintain a wormlike character similar to NaHA without bound dendrimers, those of [NH(2)]/[COO(-)] above 10 behave like rigid rods. It is concluded that the hydrogen-bonding interaction, besides the electrostatic interaction, should play an important role in the formation of the NaHA-dendrimer complexes.     

Liquid‐like structure of polymer solutions near the overlap concentration

The light scattering structure factor S( q , c) has been measured for a series of concentrations near the overlap value c* for solutions of high molecular weight poly(α‐methyl styrene) in the good solvent toluene. Scattering functions near and above overlap are characterized by a maximum as a function of scattering vector q . Scattering functions have also been calculated for these conditions using the measured second virial coefficient and radius of gyration, as reported previously for dilute solutions. The scattering function is factored into an intramolecular part that is described by a Debye function with no adjustable parameters and an intermolecular part that depends on the coil–coil pair correlation function, as suggested by Flory and Bueche. The pair correlation function is calculated using the Percus–Yevick theory of liquids and the Flory–Krigbaum potential for coil–coil interactions, as suggested by Frank Stillinger. Good agreement is obtained for the most concentrated dilute solutions, but as the overlap concentration is approached significant discrepancies are observed. The thermodynamic value of the scattering function, S(0, c), is overestimated by the theory. This discrepancy is discussed in terms of the importance of three‐body interactions, the failure of the Flory–Krigbaum potential in semidilute solutions and the limited precision of the standard protocol for calculating the measured scattering function in nondilute solutions. The observed maximum in the scattering function near overlap is not quantitatively reproduced by the theory. This discrepancy is discussed in terms of the failure of the shape of the Flory–Krigbaum potential to accurately reflect the energy of overlap for chains separated by distances near twice the radius of gyration. The mean‐field nature of the potential ignores the increased probability of interactions of linear neighboring segments. Well into the overlap region, the calculated scattering function poorly describes the observed results. The failure of the Flory–Bueche approximation in semidilute solutions is discussed as well as the effect of a changing radius of gyration as a function of concentration on the intramolecular scattering function. © 2006 Wiley Periodicals, Inc. J Polym Sci Part B: Polym Phys 44: 703–710, 2006     

Hydrodynamic interaction in polymer solutions simulated with dissipative particle dynamics

The authors analyzed extensively the dynamics of polymer chains in solutions simulated with dissipative particle dynamics (DPD), with a special focus on the potential influence of a low Schmidt number of a typical DPD fluid on the simulated polymer dynamics. It has been argued that a low Schmidt number in a DPD fluid can lead to underdevelopment of the hydrodynamic interaction in polymer solutions. The authors' analyses reveal that equilibrium polymer dynamics in dilute solution, under typical DPD simulation conditions, obey the Zimm [J. Chem. Phys. 24, 269 (1956)] model very well. With a further reduction in the Schmidt number, a deviation from the Zimm model to the Rouse model is observed. This implies that the hydrodynamic interaction between monomers is reasonably developed under typical conditions of a DPD simulation. Only when the Schmidt number is further reduced, the hydrodynamic interaction within the chains becomes underdeveloped. The screening of the hydrodynamic interaction and the excluded volume interaction as the polymer volume fraction is increased are well reproduced by the DPD simulations. The use of soft interaction between polymer beads and a low Schmidt number do not produce noticeable problems for the simulated dynamics at high concentrations, except for the entanglement effect which is not captured in the simulations.     

Hydrodynamic radius of polyethylene glycol in solution obtained by dynamic light scattering

In order to compare the size characterizations in poly(ethylene glycol) (PEG) obtained by dynamic light scattering (DLS) and small angle neutron scattering (SANS), DLS experiments were performed in various PEG solutions to ascertain the hydrodynamic radius. Data from the experiments were analyzed by using a method to eliminate effects of PEG aggregation on dynamic correlation functions. The results of the analysis were then compared to the radii of gyration reported from SANS experiments. The relation between the hydrodynamic radius, obtained by DLS, and the radius of gyration, obtained by SANS, in PEG in solution was found to be in agreement with a previously obtained relation for PEG, where the radius of gyration was found by static light scattering.     

Cutoff radius effect of the isotropic periodic sum and Wolf method in liquid-vapor interfaces of water

As a more economical but similarly accurate computation method than the Ewald sum, the isotropic periodic sum (IPS) method for nonpolar molecules (IPSn) and polar molecules (IPSp), along with the Wolf method are of interest, but the cutoff radius dependence is an important issue. To evaluate the cutoff radius effect of the three methods, a water-vapor interfacial system has been studied by molecular dynamics. The Wolf method can produce adequate results for surface tension compared to that of the Ewald sum (within 2.9%) at a long enough cutoff radius, r(c). However, the estimation of the electrostatic potential profile and dipole orientational function is poor. The Wolf method cannot estimate electrostatic configuration at r(c) ≤ L(z)∕2 (L(z) is the longest lattice of the system). We have found that the convergence of the surface tension and the electrostatic configuration of the IPSn method is faster than that of the IPSp method. Moreover, the IPSn method is most accurate among the three methods for the same cutoff radius. Furthermore, the behavior of the surface tension against the cutoff radius shows a greater difference for the IPSn and IPSp method. The surface tension of the IPSp method fluctuates and presents a similar result to that of the Ewald sum, but the surface tension for the IPSn method greatly deviates near r(c) = L(z)∕3. The cause of this deviation is the difference between the interfacial configuration of the water surface and the cutoff treatment of the IPS method. The deviation becomes insignificant far from r(c) = L(z)∕3. In spite of this shortcoming, the IPSn method gives the most accurate result in estimating the surface tension at r(c) = L(z)∕2. From all the results in this work, the IPSn and IPSp method have been found to be more accurate than the Wolf method. In conclusion, the surface tension and structure of water-vapor interface can be calculated by the IPSn method when r(c) is greater than or equal to the longest lattice of the system. The IPSp method and the Wolf method require a longer cutoff radius than the longest lattice of the system to estimate interfacial properties.     

Structural investigations of poly(amido amine) dendrimers in methanol using molecular dynamics

The molecular structure of poly(amido amine) dendrimers is investigated with molecular dynamics simulations using the Amber 7 modeling package. A method for defining residues for complex molecules is developed, and it enables the study of the effects of protonation of the primary and tertiary amines. The effects of implicit solvents versus explicit solvents as well as the pH of the solution on the molecular structure are calculated. Good agreement with experimental results for the radius of gyration measured in methanol by X‐ray scattering is observed for simulations with an explicit solvent and protonation of the primary amines. Calculations of the intramolecular atomistic pair correlation function show a dense core, as well as the presence of voids filled with the solvent inside the molecule. The primary amines (end groups) are shown to access the molecular interior by backfolding. Comparisons with experiments and other reported simulation results highlight the advantages of the approach developed here. © 2006 Wiley Periodicals, Inc. J Polym Sci Part B: Polym Phys 44: 3062–3077, 2006     

Computer simulations of diffusion and dynamics of short-chain polyelectrolytes

Brownian dynamics simulations are conducted to investigate the diffusional and dynamic properties of polyelectrolytes in dilute salt-free solutions. The polyelectrolyte molecule is represented by a bead-spring chain in a primitive model. The long-range hydrodynamic and Coulomb interactions are both taken into consideration through the Ewald summations for the first time. The major finding of our simulations is that the dependence of the long-time chain diffusivity on the Coulomb interaction strength is very different from that of the Kirkwood short-time diffusivity, which simply shows a trend nearly opposite to the chain size. When ignoring the hydrodynamic interaction (HI), the coupling effect between the chain and its counterions gives rise to a noticeable increase in the long-time diffusivity at intermediate electrostatic interaction strengths. However, the incorporation of HI suppresses this effect to a degree that one can no longer discern it. Moreover, the rotational relaxation is found to show a dependence opposite to that of the gyration radius relaxation.     

Amorphous silica modeled with truncated and screened Coulomb interactions: a molecular dynamics simulation study

We show that finite-range alternatives to the standard long-range pair potential for silica by van Beest et al. [Phys. Rev. Lett. 64, 1955 (1990)] might be used in molecular dynamics simulations. We study two such models that can be efficiently simulated since no Ewald summation is required. We first consider the Wolf method, where the Coulomb interactions are truncated at a cutoff distance rc such that the requirement of charge neutrality holds. Various static and dynamic quantities are computed and compared to results from simulations using Ewald summations. We find very good agreement for rc approximately 10 A. For lower values of rc, the long-range structure is affected which is accompanied by a slight acceleration of dynamic properties. In a second approach, the Coulomb interaction is replaced by an effective Yukawa interaction with two new parameters determined by a force fitting procedure. The same trend as for the Wolf method is seen. However, slightly larger cutoffs have to be used in order to obtain the same accuracy with respect to static and dynamic quantities as for the Wolf method.     

Melting temperature of ice Ih calculated from coexisting solid-liquid phases

We carried out molecular-dynamics simulations by using the two-phase coexistence method with the constant pressure, particle number, and enthalpy ensemble to compute the melting temperature of proton-disordered hexagonal ice I(h) at 1-bar pressure. Four models of water were considered, including the widely used TIP4P [W. L. Jorgensen, J. Chandrasekha, J. D. Madura, R. W. Impey, and M. L. Klein, J. Chem. Phys.79, 926 (1983)] and TIP5P [M. W. Mahoney and W. L. Jorgensen J. Chem. Phys.112, 8910 (2000)] models, as well as recently improved TIP4P and TIP5P models for use with Ewald techniques-the TIP4P-Ew [W. Horn, W. C. Swope, J. W. Pitera, J. C. Madura, T. J. Dick, G. L. Hura, and T. Head-Gordon, J. Chem. Phys.120, 9665 (2004)] and TIP5P-Ew [S. W. Rick, J. Chem. Phys.120, 6085 (2004)] models. The calculated melting temperature at 1 bar is T(m) = 229 +/- 1 K for the TIP4P and T(m) = 272.0 +/- 0.6 K for the TIP5P ice I(h), both are consistent with previous simulations based on free-energy methods. For the TIP4P-Ew and TIP5P-Ew models, the calculated melting temperature is T(m) = 257.0 +/- 1.1 K and T(m) = 253.9 +/- 1.1 K, respectively.     

Is the Ewald summation still necessary? Pairwise alternatives to the accepted standard for long-range electrostatics

We investigate pairwise electrostatic interaction methods and show that there are viable computationally efficient (O(N)) alternatives to the Ewald summation for typical modern molecular simulations. These methods are extended from the damped and cutoff-neutralized Coulombic sum originally proposed by Wolf et al. [J. Chem. Phys. 110, 8255 (1999)]. One of these, the damped shifted force method, shows a remarkable ability to reproduce the energetic and dynamic characteristics exhibited by simulations employing lattice summation techniques. Comparisons were performed with this and other pairwise methods against the smooth particle-mesh Ewald summation to see how well they reproduce the energetics and dynamics of a variety of molecular simulations.     

Theoretical aspects and computer simulations of flexible charged oligomers in salt-free solutions

The structural and thermodynamic properties of a model solution containing flexible charged oligomers and an equivalent number of counterions were studied by means of the canonical Monte Carlo simulation and integral equation theory. The oligomers were represented as freely jointed chains of charged hard spheres. In accordance with the primitive model of electrolyte solutions, the counterions were modeled as charged hard spheres and the solvent as a dielectric continuum. Simulations were performed for a set of model parameters, independently varying the chain length and concentration of the oligomers. Structural properties in the form of pair distribution functions were calculated as functions of model parameters. In addition, thermodynamic properties such as the excess energy of solution and the excess chemical potential of counterions were obtained. These properties were correlated with the conformational averages of oligomers as reflected in the end-to-end distances and radii of gyration obtained from the simulations. The relation with the experimental data for heats of dilution and for the activity coefficient is discussed. Finally, theories based on Wertheim's integral equation approach (product reactant Ornstein-Zernike approach) [J. Stat. Phys. 42, 477 (1986)] in the so-called polymer mean spherical and polymer hypernetted chain approximations were tested against the new and existing computer simulations. For the values of parameters examined in this study, the integral equation theory yields semiquantitative agreement with computer simulations.     

Density functional theory of liquid crystals and surface anchoring: hard Gaussian overlap-sphere and hard Gaussian overlap-surface potentials

This article applies the density functional theory to confined liquid crystals, comprised of ellipsoidal shaped particles interacting through the hard Gaussian overlap (HGO) potential. The extended restricted orientation model proposed by Moradi and co-workers [J. Phys.: Condens. Matter 17, 5625 (2005)] is used to study the surface anchoring. The excess free energy is calculated as a functional expansion of density around a reference homogeneous fluid. The pair direct correlation function (DCF) of a homogeneous HGO fluid is approximated, based on the optimized sum of Percus-Yevick and Roth DCF for hard spheres; the anisotropy introduced by means of the closest approach parameter, the expression proposed by Marko [Physica B 392, 242 (2007)] for DCF of HGO, and hard ellipsoids were used. In this study we extend an our previous work [Phys. Rev. E 72, 061706 (2005)] on the anchoring behavior of hard particle liquid crystal model, by studying the effect of changing the particle-substrate contact function instead of hard needle-wall potentials. We use the two particle-surface potentials: the HGO-sphere and the HGO-surface potentials. The average number density and order parameter profiles of a confined HGO fluid are obtained using the two particle-wall potentials. For bulk isotropic liquid, the results are in agreement with the Monte Carlo simulation of Barmes and Cleaver [Phys. Rev. E 71, 021705 (2005)]. Also, for the bulk nematic phase, the theory gives the correct density profile and order parameter between the walls.     

Shear‐Induced Diffusion in Dilute Suspensions of Charged Colloids

We propose a model for the nonequilibrium enhancement of colloidal self‐diffusion in an externally imposed shear flow in charged systems. The diffusion enhancement is calculated in terms of electrostatic, two‐body interactions between the particles in shear flow. In the high‐shear rate, low‐volume fraction limit in which the model is valid, we compare these calculations to the experiments of Qiu et al. [PRL 61, 2554 (1988)] and simulations of Chakrabarti et al. [PRE 50, R3326 (1994)] and find good agreement on scaling and magnitude to within experimental uncertainty of the electrostatic parameters.