Dissipative Particle Dynamics with an Effective Pair Potential from

Article pubs.acs.org/JPCB

Dissipative Particle Dynamics with an Effective Pair Potential from Integral Equation Theory of Molecular Liquids Alexander E. Kobryn,† Dragan Nikolić,†,‡ Olga Lyubimova,†,‡ Sergey Gusarov,† and Andriy Kovalenko*,†,‡ †

National Institute for Nanotechnology, National Research Council of Canada, 11421 Saskatchewan Drive, Edmonton, Alberta T6G 2M9, Canada ‡ Department of Mechanical Engineering, University of Alberta, Edmonton, Alberta T6G 2G8, Canada ABSTRACT: We present a method of DPD simulation based on a coarse-grained effective pair potential obtained from the DRISM-KH molecular theory of solvation. The theory is first used to calculate the radial distribution functions of all-atom solute monomers in all-atom solvent and then to invert them into an effective pair potential between coarse-grained beads such that their fluid without solvent accounts for molecular specificities and solvation effects in the all-atom system. Bonded interactions are sampled in relatively short MD of the all-atom system and modeled with best multi-Gaussian fit. Replacing the heuristically defined conservative force potential in DPD, the coarse-grained effective pair potential is free from the artificial restrictions on potential range and shape and on equal volume of solute and solvent blobs inherent in standard DPD. The procedure is flexible in specifying coarse-grained mapping and enormously increases computational efficiency by eliminating solvent. The method is validated on polystyrene chains of various length in toluene at finite concentrations for room and polystyrene glass transition temperature. It yields the chain elastic properties and diffusion coefficient in good agreement with experiment and all-atom MD simulations. DPD with coarse-grained effective pair potential is capable of predicting both structural and dynamic properties of polymer solutions and soft matter with high accuracy and computational efficiency.

1. INTRODUCTION The steady increase in the complexity of new industrial materials and the need to accurately predict their properties have resulted in rapid development of sophisticated methods of modeling on multiple length and time scales.1−20 A simple hierarchical approach combines different methods for different scales (quantum mechanical, microscopic, mesoscopic) interpreted independently. Information obtained at one level is transmitted to the next level without critical assessment or performance improvement. An advanced hierarchical approach works in a similar manner but with a feedback control. In these approaches, continuous (macro- and mesoscale) properties of materials stem from and can be controlled with their microscopic structure. To connect to large scales, molecular models representing microscopic structure must be coarsegrained (CG) to reduce detailed degrees of freedom while preserving relevant ones. The industrial importance of coarsegraining techniques to connect multiple length and time scales and to properly take into account the most notable mesoscale phenomena is increasingly appreciated.21−29 When implemented in software packages, such modeling techniques significantly accelerate the development of new materials satisfying modern industrial requirements.30,31 To interpret and guide experimental studies of structural and dynamic properties of complex systems and to predict and control the effect of underlying atomistic structure on © XXXX American Chemical Society

macroscopic properties of the system, any exploration of mesoscopic phenomena must span in length and time scales over atomistic and mesoscopic regimes. The appearance of coarse-graining as a computational framework attaining the essential physics across the range of system scales establishes an efficient mapping between subsystems on different scales. However, the loss of structural resolution in a coarse-graining procedure increases the inherent risk of omitting crucial physical properties of the system. Conveniently, as a wellestablished theoretical framework, statistical mechanics provides a mapping between disparate scales of macroscopic observables and underlying molecular interactions. Thus, it can be employed to develop new or improve existing coarsegraining techniques. In this contribution, we present and discuss a hierarchical procedure of bridging the gap between atomistic and macroscopic modeling through mesoscopic level of description. In particular, we develop a new methodology, which replaces the heuristically defined interaction potential in dissipative particle dynamics (DPD) with a CG effective pair potential obtained using molecular theory of solvation, and validate the method on a well-studied real system. We begin Received: April 23, 2014 Revised: August 26, 2014

A

dx.doi.org/10.1021/jp503981p | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

However, obtaining forces from theory is not as straightforward as correlation functions. In the fourth case, CG potentials are specified empirically in an analytical form with parameters adjusted to fit atomistic simulation data and/or experimental measurements for thermodynamic properties such as density and pressure, partitioning between phases, and solubility data. In particular, this approach was successfully used in the CG modeling of formation and rupture of micelles and spontaneous formation of lipid bilayers.34,35 Much as in the previous cases, the requirement to run all-atom MD simulation is an integral part of the approach. The CG potentials obtained can be successful in reproducing the properties to which they were fit; however, the main weakness of empirical coarse-graining techniques is that one cannot be sure if they are equally good in reproducing other properties. Furthermore, these potentials are functionals of the thermodynamic state of the system, the feature common to approaches (ii)−(iv). Therefore, these methods are sometimes regarded as nonrigorous. Alternatively, CG modeling can be dynamics based, which means that evolution of the system is described with some equation for dynamic properties. This is typically done through the Langevin equation45 in which the friction coefficient is tuned until the CG dynamic properties match the atomistic ones. Although the formalism can be applied to a wide variety of systems, its major limitation is the inability to reproduce structural properties. Finally, coarse-graining in time describes evolution over long time scales while preserving description on microscopic level. This is achieved by developing new or combining existing algorithms allowing extension of accessible time scales as compared to standard all-atom simulations. Because CG modeling aims at quickly yielding properties of a large microscopic system on a very long time scale, a theoretical approach is highly desirable. The requirement to begin modeling with running MC or MD simulation dramatically increases computational cost. Among the existing theories dealing with equilibrium correlation functions, of practical interest for coarse-graining purposes is the Ornstein−Zernike (OZ)-type integral equation theory of liquids and solutions.46,47 The idea of using the one-component OZ equation complemented with the hypernetted chain (HNC) closure to construct the effective potential between CG polymers has been advocated in ref 48. In that study, to enable modeling of large-scale mixtures of polymers and colloids, it was suggested that each polymer is replaced by a single particle interacting with other polymer particle via an effective pair potential based on a center-of-mass (CM), midpoint, or end-point representation. A comparison to the reference microscopic self-avoiding walks (SAW) system directly simulated by employing the MC pivot algorithm together with translational moves and configurational bias MC algorithms showed a good agreement and a surprisingly weak density dependence of the effective potential. The latter observation can be understood from the effect of density-independent many-body interactions, which for low and moderate densities are dominated by three-body interactions.49 The molecular version of OZ theory in the interaction site formalism, called reference interaction site model (RISM),50 relates site−site RDFs to direct correlation functions (DCFs) having the asymptotics of site−site interaction potentials. A proof of concept that a coarse-graining model reduction strategy based on integral equation theory is feasible within just the RISM-HNC approximation was given in ref 51 where it

with outlining the existing coarse-graining approaches and pointing out their strong and weak sides. The current state of multiscale modeling based on coarsegraining can be conventionally classified into the following categories:32−35 (i) rigorous coarse-graining methods; (ii) coarse-graining by matching correlation functions; (iii) coarse-graining by matching forces; and (iv) empirical coarsegraining. In the first case, the partition function of the fully atomistic system is rigorously mapped onto that of the CG system, which ensures that the new system accurately reproduces all of the equilibrium thermodynamic properties of the initial one. This is the main advantage of approach (i). It is achieved through construction of the effective Hamiltonian for the CG system in which part of the interaction terms is replaced with an effective interaction.36,37 The method turned out to be useful to model phase behavior of some neutral asymmetric mixtures, but is not applicable to systems with long-range interactions or with structured molecules, mostly because calculation of the effective interaction term becomes computationally expensive or impractical. In the second case, it is customary to characterize CG species with a potential of mean force (PMF) defined as the Boltzmann inversion of a two-body radial distribution function (RDF) between interaction sites in the fully atomistic system. The CG system is supposed to reproduce a set of correlation functions of the initial fully atomistic system, for example, a set of RDFs. The matching procedure assumes that correlation functions of the initial system can be obtained from either Monte Carlo (MC) molecular simulation, molecular dynamics (MD) simulation, or analytically from theory. This method is very popular to deal with long-range electrostatic interactions or complex intramolecular degrees of freedom. It works equally well for systems with short-range interactions. Note that inversion of an RDF in a dense fluid into a CG interaction potential between two CG spherical particles is unique. However, it does not guarantee that all of the equilibrium thermodynamic properties of the atomistic system are reproduced by the CG one,38 as coarse-graining contracts molecular degrees of freedom. (The disparity can be minimized by optimizing the mapping of all-atom species into CG ones, which is not unique and is usually chosen for geometric considerations and convenience.) Furthermore, if a thermodynamically consistent model is required, a possibility of mixed CG strategies that try to match simultaneously not only the pair distribution function but also some intramolecular three-body and four-body distributions and some thermodynamic property, for instance, by constraining the virial pressure to be equal to the pressure of the microscopic model, is discussed.39,40 We emphasize that the additional constraint on pressure consistency makes the matching problem overcomplete and so changes it into an optimization problem. In the third case, the CG potential is determined to fit the forces on the mapping interaction sites of the CG model to those obtained from all-atom MD simulation.41−44 An additional condition to be satisfied in such a procedure is that the virial of the CG system matches that of the atomistic simulation. Similarly to approach (ii), the effective potentials so obtained are temperature and density dependent, but can help better predict thermodynamic properties such as pressure or density. It is also well transferable to different system sizes. In general, this approach is applicable to arbitrary systems for which all-atom MD simulations are possible in principle. B



Article

tration of the all-atom system. The resulting CG system without solvent and/or counterions is a fluid of CG species at a considerably reduced packing fraction, interacting with an essentially short-range effective potential, which makes it well suitable for DPD. Note that the method is highly flexible in a choice of species for structural coarse-graining and for contraction in a multicomponent system. This provides a dramatic speedup and yet high accuracy of the modeling with focus on a particular component in effective potential DPD simulation of a CG contracted system. To this end, we propose a new methodology that matches the RDFs of a fluid of structurally CG species to the site−site RDFs of the full multicomponent system of atomic species by using the DRISM-KH integral equation theory of molecular liquids. This Article is organized as follows. In section 2, we describe how the so-called conservative force of DPD is replaced in our new methodology with an effective force comprising bonded CG interactions sampled and modeled separately and a nonbonded CG effective pair potential obtained from the DRISM-KH molecular theory of solvation, and how the molecular effect of contracted solvent (and/or counterions) is treated in this approach without loss of modeling accuracy. For convenience, the standard DPD approach along with its limitations such as an empirical definition of the so-called conservative force are briefly reviewed in Appendix A. The treatment of bonded interactions in a CG polymer chain is described in Appendix B. To support the derivation of the integral equations for a nonbonded CG effective pair potential, the main equations of the DRISM-KH molecular theory of solvation are shown in Appendix C. The flowchart diagram for DPD with CG effective pair potential from DRISM-KH is put into Appendix D. In section 3, we provide the CG model of polystyrene and the setup of the DRISM-KH calculation and of the MD and DPD simulations to validate the method on CG modeling of polystyrene in toluene solution at different temperatures and compositions, including different polystyrene tacticities. Results and discussion of this validation case are presented in section 4, including the effective pair potential, the bonded CG interaction potentials, and the DPD results for the polymer stiffness and the diffusion coefficient. The calculated polymer stiffness is compared to our all-atom MD simulation data, and the results for the diffusion coefficient are benchmarked against experimental data. In section 5, we summarize the work and outline future studies.

was tested on a simple model system of rigid cyclohexanol as an example of an amphiphilic solute system containing both hydrophobic and hydrophilic regions, solvated in TIP3P water under ambient conditions. To yield an effective potential between a reduced set of atomic sites and the solvent, the RISM-HNC equation was successfully inverted in a numerically stable way, thus showing that a bottom-up approach to reduced model development starting from established molecular force fields is possible without referring to iterative simulation-based compression.51 In the present work, to derive an effective pair potential suitable for DPD simulations that coarse-grains selected species (e.g., polymer beads) and contracts all other solution components (e.g., solvent and counterions), we employ the dielectrically consistent (DRISM) theory52,53 complemented with the Kovalenko−Hirata (KH) closure relation.54−56 The DRISM-KH theory theory is analytically renormalized to ensure proper electrostatic asymptotics of site−site RDFs in electrolyte solution at a finite concentration and is capable of treating various liquids, mixtures, and solutions comprising nonpolar, polar, and ionic molecular species of a given composition in a wide range of thermodynamic conditions and local environment of nanoporous materials, soft matter, and biomolecular systems.50,55−59 In particular, it provided an insight into the experimentally observed structural transitions and related thermodynamic anomalies in water−alcohol mixtures60,61 and into the microscopic structure of oligomeric polyelectrolyte gelators in different solvents.57 It is free from the limitations inherent in heuristic theories used in standard DPD simulations, including the restrictions on the potential range and shape and the requirement of equal volume of solute and solvent CG beads inherent in Flory−Huggins theory. In terms of the coarse-graining approaches outlined above, our method is spatial structure based and thus falls into category (ii). As distinct, it employs integral equation theory of molecular liquids for both the fully atomistic system and the CG one, with the advantages of the analytical treatment of long-range asymptotics of correlation functions and the absence of noise inherent in molecular simulations. The method is also much faster than molecular simulation of all-atom RDFs and iterative Boltzmann inversion based on a sequence of molecular simulations, which becomes critical for systems with strongly associating species, hydrogen bonding, and slow exchange of solvent and ions, such as oligomeric polyelectrolyte gelator with different counterions in different solvents.57 First, the DRISM-KH theory is applied to obtain the site− site RDFs for the full all-atom system. Next, the selected RDFs are mapped into those of a fluid comprising just the CG species at the same partial number densities, with all of the other components contracted. Next, the DRISM-KH equations for the fluid of CG species are used to invert the RDFs into effective pair potentials between CG beads. If the atomic species selected for coarse-graining are ionic, their counterions can be contracted together with solvent and other components of the solution. In this case, all of the electrostatic interactions of the all-atom system are mapped onto the CG effective pair potential, which acquires a screened electrostatic asymptotics with the dielectric screening by all polar species in the solution and the Debye screening by all ionic species except for those being coarse-grained. Note that the CG effective pair potential with this “partial” screening Yukawa asymptotics inserted in the DRISM-KH equations leads to the RDFs of the CG fluid reproducing the Debye screening by the full ionic concen-

2. DPD WITH AN EFFECTIVE PAIR POTENTIAL DPD is a stochastic dynamics method originally developed for simulations of hydrodynamic phenomena in simple and complex fluids62 in which individual mesoscale particles (beads) represent a set of atoms or regions of fluids. Because the original formulation did not satisfy the fluctuation− dissipation theorem, a modified DPD algorithm63,64 was formulated to be fully consistent with statistical mechanics, thereby establishing DPD as a valid method to study the dynamics of mesoscopic systems. It has been applied previously to various problems in rheology,65−67 material science,68−71 biochemistry,72,73 molecular biology to model membranes,74−76 vesicles,77−79 and micellar systems,80−82 etc. A detailed formulation of DPD widely used at present can be found, for example, in refs 83−86. For convenience and comparison with our modification, we present the DPD basics in Appendix A and also outline the treatment of bonded CG interaction potentials in DPD in Appendix B. C



Article

2.1. Effective Solvent Treatment in DPD. Applications of DPD include soft matter structures of large assemblies of macromolecules with various cofactors in solution and fluid flow in complex geometries. A substantial part of computing time, especially for dilute solutions, is spent on describing motion and exchange statistics of solvent molecules. Regardless of interparticle potentials in the whole system comprising solute and solvent, it is perfectly legitimate to use a DPD thermostat that treats solvent implicitly as an effective continuous medium.87 Such a treatment dramatically decreases the computing time and yet retains the main benefit of the DPD technique: correctly reproducing the hydrodynamic correlations while conserving the system momentum. To extend this idea to realistic systems, it is essential to derive a suitable interaction potential between DPD beads that not only accounts for molecular packing effects88,89 but also represents chemical functionalities of the polymer and solvent. To develop such an effective pair potential has been another major motivation for the present study. 2.2. Development of an Effective Pair Potential. Effective interaction potentials of species in solution can be deduced by analysis of RDFs obtained from either experiment or molecular simulations. In particular, the problem of reconstruction of effective interaction potentials from RDFs was addressed by using the inverse Monte Carlo method combined with the CG implicit solvent model90 to study ion− DNA interactions in solution.91 This iterative approach relies on a statistical-mechanical argument that a local correction to effective interaction potentials propagates in RDFs to all distances and that the iterative difference in ionic interaction potentials is solvent mediated.92 In this work, we propose to describe macromolecules in solution with an approach of effective potential DPD in which an effective pair potential obtained from molecular theory of solvation accounts for both a CG interaction potential and effective solvation forces between solute macromolecules, statistically mechanically averaged over arrangements of allatom solvent molecules in the solvation shells. Thus, a CG polymer system described with such an effective pair potential contains essential information on molecular geometry, chemical specificities, and solvation structure of both polymer and solvent. In our effective potential DPD approach, the conservative force FCij of the conventional DPD force field in eq A1 is replaced with a pairwise-additive conservative force Fijeff corresponding to the effective pair potential between nonbonded CG beads of solute macromolecules ϕeff nb that comprises both the interaction potentials and the solvent effect. The momentum-conserving dissipative force FDij and random force FRij are kept to constitute a thermostat, and so the resulting force becomes Fi =

∑ (Fijeff j≠i

Specifically, the site−site RDFs for the all-atom model of solute monomers in all-atom solvent are obtained from the integral equations of the dielectrically consistent reference interaction site model (DRISM)52,53 complemented with the Kovalenko−Hirata (KH) closure relation,54−56 which are shown in Appendix C. The KH closure approximation consistently accounts for both electrostatic and nonpolar features of solvation (associative and steric effects) in simple and complex liquids, ionic liquids, and nonelectrolyte and electrolyte solutions in various chemical systems.55,56 In particular, the DRISM-KH approach provided an insight into the structural transitions and related thermodynamic anomalies in the formation of micromicelles and the transition from the tetrahedral to zigzag hydrogen-bonding network in water− alcohol mixtures in the whole range of concentrations.60,61 As another example, DRISM-KH was an essential part of a multiscale approach to the microscopic structure of gels formed by oligomeric polyelectrolyte gelators in different solvents.57 Converging the DRISM-KH integral equations by using the modified direct inversion in iterated space (MDIIS) accelerated numerical solver55,93,94 yields the site−site RDFs gαγ(r) for the all-atom solution of monomers. We next require that the CG bead−bead RDFs g̃αγ(r) of a fluid of CG species at the same CG site number density ρ̃α = ρα interacting with the nonbonded CG effective pair potential φ̃ αγ(r) exactly match the all-atom RDFs between interaction sites used in coarse-graining mapping, g̃αγ(r) = gαγ(r) for α, γ ∈ CG. The DRISM-KH integral eqs C1 and C10 are now used for the fluid of CG species with the corresponding CG total and direct correlation functions h̃αγ(r) and c̃αγ(r), where all site indices now enumerate interaction sites of only CG species (more than one in a general case). The CG intramolecular ̃ (r) matrix ω̃ αγ(r) and the renormalized dielectric correction ζαγ have the same form C4 and C5 but with the site separations and the parameters in eqs C6−C9 determined by the geometry and site charges of CG species. The DRISM-KH inversion of the CG RDFs g̃αγ(r) into the effective pair potential φ̃ αγ(r) is readily performed in two single steps. The CG direct correlation function c̃αγ(r) is calculated from the DRISM equation in the reciprocal space as ̃ (k) − ζμν ̃ (k)ρ ̃ ]−1 *[hμν ̃ (k)] cαγ ̃ (k) + ζαμ ̃ (k) = [ωαμ μ ̃ (k)]−1 *[ω̃νγ(k) + ρν̃ hνγ

(2)

where the asterisk * denotes here summation over repeating CG site indices, h̃αγ(r) = g̃αγ(r) − 1, and transformation between the direct and reciprocal space is performed with the analytical treatment of the electrostatic asymptotics55,56 and discrete sine transform for the remaining short-range terms. The effective pair potential φ̃ αγ(r) is then calculated from the KH closure C10 as

+ FijD + FijR )riĵ (1)

φ̃αγ(r ) = kBT[gαγ ̃ (r ) − 1 − ln gαγ ̃ (r ) − cαγ ̃ (r )]

Extracting long-range asymptotics of radial distribution functions from molecular simulation data and filtering out noise inherent in molecular simulations pose a challenge for inverse methods of potential reconstruction from distribution functions, including iterative Boltzmann inversion. To avoid these issues, we employ statistical-mechanical, molecular theory of solvation,50 first to obtain the distribution functions of the all-atom solution model and then to invert them into the effective pair potential of the CG fluid.

for gαγ ̃ (r ) ≤ 1 = −kBTcαγ ̃ (r ) for gαγ ̃ (r ) > 1

(3)

where k B is the Boltzmann constant and T is the thermodynamic temperature. For all-atom solution with ionic species, the Coulomb interaction potentials can be contracted to Yukawa-type CG D



Article

solutions of atactic polystyrene in toluene. Atactic polystyrene is chosen for this study because it is a commodity polymer with a broad range of end-uses and one of the most studied both experimentally and numerically.95−98 Figure 1 illustrates the

effective pair potentials in the CG system. If the atomic species selected for coarse-graining are ionic, their counterions can be contracted together with solvent and other components of the solution. (For a molten salt, counterions can be so contracted together with coarse-graining of molecular or polyelectrolyte ions.) In this case, all of the electrostatic interactions of the allatom system are mapped onto the CG effective pair potential, which acquires a screened electrostatic asymptotics: eff 2 ϕnb (r ) ≈ qCG

exp( −r /lDeff ) εr

(4)

with a “partial” Debye length: lDeff

⎛ ⎞−1/2 2 ⎜ = ⎜4π ∑ qa ρa /(εkBT )⎟⎟ ⎝ a ≠ CG ⎠

(5)

of the ionic screening by all species a with ionic charge qa at concentration ρa except for the ionic species being coarsegrained. It can be readily shown that the CG effective pair potential with this “partial” screening Yukawa asymptotics inserted in the DRISM-KH equations for the CG fluid produces the RDF and the corresponding potential of mean force (PMF) possessing the electrostatic asymptotics: wnb ̃ (r ) ≈ 1 −

2 qCG exp( −r /lD) kBT εr

(6) Figure 1. Atomistic structures of the polystyrene stereoisomers: isotactic (a), syndiotactic (b), and atactic (c). Red “●” indicate carbon atoms C(3) used for mapping between the atomistic and CG models of polystyrene. Coarse-grained models (d): isotactic polystyrene bead (iPS) centered at interaction site C(3), and toluene bead with interaction site C(9). Encircled areas mark the iPS bead centered at interaction site C(3) and the toluene bead with interaction site C(9).

with the Debye length lD due to the full ionic concentration of the all-atom system, including the ionic solute. Furthermore, for a single-site CG bead (which is the validation case we present below), the DRISM-KH eqs C1 and C10 simplify to the OZ-KH ones. The effective pair potential between nonbonded CG beads ϕeff nb(r) is then obtained by the OZ-KH inversion of the bead−bead RDF g̃(r) of the CG fluid at density ρ̃: c (̃ k) = h(̃ k)/[1 + ρ ̃h(̃ k)]

scheme of mapping of the all-atom structure of polystyrene into the corresponding CG beads. It is similar to that proposed by Sun and Faller,99 with the difference that our mapping site coincides not with the center of mass of the polystyrene repeat unit but with the nearby C(3) atomic site. Statistical correlations of internal degrees of freedom depend on a choice of mapping sites such are C(2) carbon atoms in Figure 1d, and further information on possible mapping strategies can be found in the review article by Karimi-Varzaneh et al.98 The 3D solvation structure of an atactic styrene trimer (solute) in toluene solvent is visualized in Figure 2 with the 3D isosurface of the PMF of toluene solvent site C(9) around the solute at infinite dilution in ambient toluene. The potential isovalue corresponds to the thermal energy at room temperature. Note that as distinct from the conventional DPD approach, our method needs no coarse-graining of solvent, as the effective pair potential between polymer beads (styrene−styrene) already accounts for all-atom solvent-mediated interactions (toluenemediated contribution). 3.2. MD Simulation Setup. Dense environment significantly affects flexible intramolecular bonds, with a nontrivial difference between the effect of solvents of different molecular nature, as observed, for instance, for oligomeric electrolyte gelators in a set of solvents.57 To sample the bonded CG interaction potentials of stretching, bending, and torsion in the present system, we performed all-atom MD simulations of a single atactic polystyrene chain of 60 repeat units solvated in toluene. Note that for better thermodynamic consistency, the

(7)

eff ϕnb (r ) = kBT[g (̃ r ) − 1 − ln g (̃ r ) − c (̃ r )] for g (̃ r ) ≤ 1

= −kBTc (̃ r ) for g (̃ r ) > 1 (8)

Note that the inversion of the CG RDFs through eqs 2,3 or 7,8 poses no numerical problems, as all of the correlation functions are obtained in the framework of the DRISM-KH theory with the asymptotics specified analytically. The effective pair potential ϕeff nb(r) so obtained contains the bead−bead direct interaction potential, including the bead−bead core repulsion appearing in the form 8 through the exponential decay of the all-atom RDF gaa(r) in the repulsive core region. ϕeff nb(r) also represents solvation effects appearing, in particular, as shortrange oscillations. gaa(r) is also expected to exhibit asymptotic decay at large bead−bead separation. The effective pair eff potential ϕeff nb(r) and its force F (r) are tabulated on a uniform radial grid and then used as an input in DPD simulations. The method to introduce and calculate the effective pair potential ϕeff nb(r) is the central result of this Article.

3. CALCULATION SETUP FOR VALIDATION ON POLYSTYRENE IN TOLUENE SOLUTION 3.1. CG Model of Polystyrene. To validate the proposed methodology, we carried out effective potential DPD for dilute E



Article

NPT ensemble MD simulation runs. The dielectric constant of the mixture of styrene monomers and toluene was approximated using Looyenga’s phenomenological expression.111 The DRISM-KH integral eqs C1−C10 were discretized on a uniform radial grid of 8192 nodes with 0.05 Å resolution and converged to a relative root-mean-square accuracy of 10−8 by using the modified direct inversion in iterated space (MDIIS) accelerated numerical solver.55,93,94 3.4. DPD Simulation Setup. We carried out mesoscopic simulations using the large-scale atomic/molecular massively parallel simulator (LAMMPS)112 for a single CG polystyrene chain in a nonperiodic cubic box. Depending on the polystyrene chain length, the box size ranges from 34 to 252 nm, and is selected so as to keep the overall bead density of about 1.25 × 10−4 nm−3. Manipulation with topology information for polystyrene chains and automated generation of input files for LAMMPS was done with the TopoTools VMD plugin.113 Initial equilibration of each CG polystyrene chain starts with generating a Gaussian ensemble of velocities with a variance scaled to the glass transition temperature such that both the linear and the angular momentum are zeroed. We then performed iterative relaxation of the system with the stopping tolerances of 10−4 and 10−6 for the energy and force evaluations, respectively. The minimization target function comprises the energy of the nonbonded effective and bonded interactions (Figures 3 and 4, respectively), such that all

Figure 2. 3D map of the potential of mean force (PMF) between molecules between the styrene trimer immersed in ambient toluene and a toluene molecule with orientationally averaged around interaction site C(9), obtained from the 3D-RISM-KH molecular theory of solvation. Isosurface (in green) of PMF = 0.59 kcal/mol (1 kT). The 3D map is calculated using the 3D-RISM-Dock algorithm implemented in the AutoDock package.141 The molecular structure and atomic site nomenclature are given in Figure 1d.

sampling points over thermodynamic states have to adequately represent the system phase states in the range of interest. In the present study, we limit ourselves to two specific conditions: at ambient pressure and temperature, as well as the polystyrene glass transition temperature.100 The simulations were carried out in the NPT ensemble by using Version 4.5.4 of the GROMACS molecular dynamics package.101−105 Temperature coupling was done by velocity rescaling106 with a stochastic term τt = 0.1 ps, and pressure regulation107 was achieved by using the Parrinello−Rahman barostat108 with τt = 2 ps. Polystyrene and toluene molecules were modeled with the OPLS-AA force field109,110 as implemented in the GROMACS package. A cutoff of 14 Å without a switch function or longrange dispersion correction was applied in the calculation of nonbonded interactions. The polystyrene chain of 60 beads was placed in the simulation box of size 60 × 60 × 120 Å filled with 2583 toluene molecules. The average physical density of the solution of polystyrene in toluene calculated over the 10 ns production run was 862.282 kg/m3 at T = 298 K and 766.884 kg/m3 at Tg = 380 K. Each simulation lasted for 20 ns (10 ns equilibration and 10 ns production run) with 2 fs time steps and neighbor list updates every five steps. At that point, we considered that the polystyrene configurations sample the intrachain conformational space well enough, and used them further as an input in mesoscopic simulations with a DPD thermostat. To validate the results of effective potential DPD simulations, in particular, for the stiffness of a solvated polystyrene chain, we also performed three MD simulations for the syndiotactic, atactic, and isotactic polystyrene stereoisomers of a single polystyrene chain of 25 beads solvated in ambient toluene, with the MD simulation setup the same as above. 3.3. DRISM-KH Setup. The molecular geometries, partial charges, and Lennard-Jones potential parameters of interaction sites in the all-atom models of styrene and ambient toluene for the DRISM-KH calculations were taken from the OPLS-AA force field.109,110 The calculations were performed for toluene solutions of styrene monomers at 0.48 wt % concentration, both at room temperature and at the polystyrene glass transition temperature. The number densities of styrene monomers and toluene were taken from the results of the

Figure 3. Results of the DRISM-KH theory for 0.48 wt % solution of styrene monomers in toluene at room temperature. Part (a): Radial distribution functions between toluene interaction site C(9) and styrene interaction sites C(2) and C(3) typically used for coarsegraining (green and blue curves, respectively), and between styrene interaction sites C(3) (red curve). Part (b): Potential of mean force w̃ nb(r) of CG styrene species (red curve) fitted to that of interaction sites C(3) of the all-atom styrene model, obtained from DRISM-KH. Effective pair potential of the nonbonded interaction of CG styrene beads ϕeff nb(r) (cyan curve).

improper interactions and external degrees of freedom are set to zero. In the first 2 ns of mesoscopic simulations, we performed a series of energy minimizations followed by trajectory relaxation in the microcanonical ensemble and scaling of the bead velocities using the Langevin thermostat. Updates of the bead positions and velocities and of the nearestF



Article

4. RESULTS AND DISCUSSION 4.1. Effective Pair Potential. We first applied the DRISMKH theory to the all-atom model of 0.48 wt % styrene monomers in ambient toluene solution, thus eliminating the conformation specific bonded interactions present in polystyrene chains and focusing only on the nonbonded component of the effective interactions. Figure 3a shows the styrene− styrene (solute−solute, or “uu”) RDF guu αγ(r) and the styrene− toluene (solute−solvent, or “uv”) RDFs guu αγ(r) for the toluene solution of styrene monomers. We then set the RDF of the fluid of CG styrene monomers g̃nb(r) equal to the RDF between sites C(3) of styrene monomers in toluene solution, and inverted it in two noniterative steps through the OZ-KH integral eqs 7,8 into the effective pair potential ϕeff nb(r) between CG monomers. Figure 3b depicts the effective pair potential ϕeff nb(r) between the CG styrene monomers in comparison with the potential of mean force obtained by the Boltzmann inversion of the CG fluid RDF as w̃ nb(r) = ln[g̃nb(r)]. The shape of the effective pair potential in general follows that of the PMF, including the positions of the first and second minima and the barrier. As compared to the effective interaction potential, the PMF has a somewhat larger oscillation amplitude and is steeper in the repulsive core region, which is a multibody effect in the fluid of CG monomers. 4.2. Results of the MD Simulations. For the two all-atom MD simulations of the atactic polystyrene chain of 60 repeat units in toluene at temperature T = 298 and 380 K, Figure 4 presents the bond stretching, bending, and torsion intramolecular properties characterized with the probability distributions Pstr(S ) of bond length S [part (a)], Pbend(θ) of bond angle θ [part (c)], and Ptors(φ) of dihedral angle φ [part (e)]. (The intramolecular variables are specified between polystyrene sites C(3) used for coarse-graining.) We validated the assumption of uncorrelated intramolecular degrees of freedom (B1) by checking that the CG model reproduces each of the intramolecular distributions obtained in atomistic sampling. Accordingly, the bonded CG interaction potential ϕeff bond is decomposed, following eq B4, into the stretching, bending, and torsion PMF terms dependent on bond length S , bending angle θ, and dihedral angle φ. Figure 4b, d, and f shows the PMF terms ϕstr(S ), ϕbend(θ), and ϕtors(φ) obtained by Boltzmann inversion of the smoothed probability distribution functions, eqs B3a−B3c. Transferability of the intramolecular PMF obtained by single-temperature sampling of the intrachain probability distribution functions to different temperatures can be achieved either by applying analytical temperature-dependent scaling factors114 or by interpolating over a set of reference temperatures.115 Note that only the torsional distribution and the corresponding PMF term exhibit considerable variations with temperature, bringing about local minima (Figure 4e and f), while the temperature variations in the bond length and angle can be neglected (Figure 4a−d). 4.3. Elastic Properties. To validate the intrinsic stiffness of a single polystyrene chain along its backbone, we calculated the mean square distance ⟨R2(N)⟩ between two arbitrary units separated by N-unit long chain segments. Trajectories in atomistic simulations are analyzed for C(3) sites, and the results obtained are illustrated in Figure 5 for 25-unit long polystyrene chains with the three different tacticities. Figure 5 also shows the ⟨R2(N)⟩ distances obtained by analyzing the trajectories from CG simulations with C(3) sites being mapping points (centers of beads). For each tacticity, we

Figure 4. Probability distributions of the bond length (a), bond angle (c), and dihedral angle (e), as well as the corresponding Boltzmanninverted PMFs (b, d, f), in an atactic polystyrene chain of 60 repeat units. Best fit with multi-Gaussian expansion (solid curves) to all-atom MD simulations (○) of the polystyrene chain in toluene solution at room temperature (blue curves and symbols) and polystyrene glass transition temperature (red curves and symbols).

neighbor lists were made every 10 fs. We then followed the DPD production run of 10 μs with the 10 fs integration step. G



Article

Figure 6. Mass dependence of the diffusion coefficient for single CG atactic polystyrene chain in toluene. Experimental measurements: ●116 and ○.120 Gray dashed line shows best fit of the experimental diffusion coefficient D0 to the power law dependence on the molar mass Mw, which works remarkably well for medium to large Mw but somewhat deviates for short chains due to chemically different end segments.119 Present study with the effective pair potential: averages of 10 and four independent DPD runs (red and blue “+”, respectively). The simulation results follow the fitting line at all Mw.

Figure 5. Stiffness of a single syndiotactic, atactic, and isotactic polystyrene chain of 25 beads solvated in ambient toluene. Trajectories of the effective potential DPD and all-atom MD simulations are analyzed with respect to the C(3) mapping sites shown in Figure 1d.

generated 10 initial configurations of the CG systems by mapping the equilibrated atomistic structures, and further subjected these to 10 μs DPD simulations with the 10 fs integration step. We sampled CG chain configurations by taking snapshots every 100 ns to characterize intrachain distances with the mean square separation ⟨R2(N)⟩ between the centers of any two beads connected via N + 1 links. As is evident from Figure 5, the proposed CG model with an effective pair potential accurately reproduces the intrinsic stiffness of the atomistic model with explicit solvent. Moreover, the observed variance in ⟨R2(N)⟩ between CG chains with different tacticity is sufficiently small. Below we focus on atactic polystyrene, as it is the only commercially important form of polystyrene whose physical behavior in dilute and semidilute solutions is studied experimentally in more detail.116 4.4. Diffusion Coefficient. The mutual diffusion coefficient of polystyrene in toluene D0 was calculated from the mean square displacement ⟨|ΔR|2⟩ = ⟨|R(t) − R(0)|2⟩ of the center of mass of each chain, averaged over all of the chains and simulation times (6D0 = limΔt→∞⟨|ΔR|2⟩Δt−1). Red and blue “+” in Figure 6 show the results for the diffusion coefficient averaged in 10 and 4 independent runs of effective potential DPD simulations, respectively. To avoid the speedup of collective diffusion for polymers in good solvent within the dilute regime, all of the effective potential DPD simulations were performed using single atactic polystyrene chains of various length. Experimental data (●) are taken from Table 8 of Rauch and Köhler,116 and the dashed straight line represents the power law fit D0 = 61(5) × Mw−0.58(2). Because the power scaling law is expected to hold only in the limit of long chains with Mw > 10 kg/mol, the experimental data for short chains of weight 0.502 and 2.41 kg/mol from Schäfer117 and 6.1 kg/mol from Köhler et al.118 are not included in the fit. There is a difference between experimental measurements and the power law fit for polymeric chains with molar mass below 10 kg/mol. Rauch and Köhler119 attribute this difference to end group effects, which can lower the diffusion coefficient of the styrene dimer in toluene by as much as 25%, as compared to longer oligomers. As seen from Figure 6, the calculated diffusion coefficient D0 as a function of Mw closely follows the power law fit at all values of Mw. The CG model we used has one type of beads for the entire polymer chain, and therefore neglects end group effects in shorter atactic polystyrene oligomers and

predicts D0 following the scaling power law up to low molar masses. The result deviates from the experimental findings of Rauch and Köhler116 and the unconfined atomistic MD result of Bayramoglu and Faller120 (shown with an ○ in Figure 6). However, we note that Bayramoglu and Faller120 report the confined value of the diffusion coefficient D0 to be 5.5 times smaller than in unconfined atomistic simulations.

5. SUMMARY We have introduced a new method of coarse-graining of soft matter systems in solution for DPD simulation. It is built on the notion of an effective pair potential of interaction between CG beads of solute macromolecules that comprises both a nonbonded CG interaction potential and effective solvation forces, obtained from integral equation theory of molecular liquids. We define the effective pair potential as an interaction potential between two nonbonded CG beads that for a fluid of CG solute macromolecules alone results in the RDF (or equivalently its Boltzmann inversion into PMF) reproducing the corresponding site−site RDF (or PMF) of the solution comprising both all-atom solute and all-atom solvent. We begin with treating the motion modes associated with intrachain bond stretching, bending, and torsion as mutually independent and much stronger in comparison to solvent-mediated and multibody effects and thus as effectively separate from the nonbonded CG effective pair potential. The bonded CG interaction potentials of stretching, bending, and torsion are sampled in a relatively short standard MD simulation of the full all-atom system, modeled with best multi-Gaussian fit to eliminate simulation noise, and used afterward. First, we calculate the site−site RDFs for the system of all-atom solute monomers in all-atom solvent from the DRISM-KH integral equation theory of molecular liquids, or molecular theory of solvation. We then map the all-atom RDF between the sites selected for coarse-graining into the RDF of a fluid of CG solute monomers at the number density corresponding to the all-atom solute component, and next calculate the effective pair potential in two single noniterative steps from the DRISM-KH H



Article

integral equations for the fluid of CG monomers. (Note that for CG beads with a single mapping site the DRISM-KH equations reduce to the OZ-KH ones.) The inversion of the RDFs into the effective pair potential in the CG fluid through the DRISMKH equations does not pose any ill-conditioning or statistical noise problems arising for RDF data from molecular simulations, as the all-atom RDFs obtained from the DRISMKH theory have both the long-range asymptotics and the core repulsion exponent specified analytically. The above analytical treatment of asymptotics of nonbonded correlations constitutes an important advantage of our method, as compared to simulation approaches to obtain CG potentials, including iterative Boltzmann inversion. This is particularly important for soft matter systems with electrostatic interactions. Moreover, our technique is suitable for fast throughput modeling with quick recalculation of the CG effective pair potential for a large set of temperatures, densities, and compositions, thus avoiding the problem of transferability inherent in CG potentials. It should be noted, however, that because coarse-graining contracts molecular degrees of freedom, this method generates a CG effective pair potential dependent on thermodynamic state, which in general is not thermodynamically self-consistent, that is, might not accurately reproduce thermodynamic properties further to the structure and dynamics of the system.38 This cannot be eliminated but can be minimized by optimizing the CG mapping, which is not unique and is usually chosen for geometric considerations and convenience. With the bonded CG interaction potentials inferred from atomistic MD sampling of local conformations, we finally use the nonbonded CG effective pair potential at the next level of multiscale modeling, in DPD simulation, which offers several notable advantages. First, a CG effective pair potential so obtained based on the first-principle foundations of the statistical-mechanical theory of liquids and solutions is free from the artificial limitations of an heuristically defined conservative force potential in the standard DPD approach based on Flory−Huggins theory. The latter typically restricts the potential range and shape and requires CG beads of solute and solvent to be of the same volume. Second, the effective pair potential includes solvation effects with full account for molecular specificities of both solute and solvent. Calculation of the CG effective pair potential by converging the DRISM-KH integral equations is a relatively fast procedure that is easily repeated to recalculate the effective pair potential for each given set of thermodynamic parameters (temperature, density, composition). This ensures that the CG model obtained with the proposed approach closely reflects chemical specificities of its parent atomistic force field and easily handles any changes in molecular structure and composition of solute and solvent. For polymeric solute, the polymer chain tacticity and length dispersion can also be easily taken into account in each particular case. Finally, a substantial computational speedup is achieved in the proposed approach of effective potential DPD because the solvent is contracted and in general the remaining fluid of CG solute is of a significantly lower packing fraction. We illustrated the method on an example of polystyrene solution in toluene, and have run effective potential DPD to describe the dynamical properties of polystyrene in dependence on its degree of polymerization. The bonded interaction potentials of CG chains were inferred from the atomistic MD sampling of local conformations of a polymer chain in toluene at infinite dilution. The stiffness of a single CG polymer chain

obtained from effective potential DPD is in a very good agreement with that from all-atom MD simulation for all of the polystyrene tacticity types and chain lengths considered. We were able to efficiently run simulations for chain length ranging from 5 to 2001 units and to observe long time scale phenomena such as diffusion in dilute solutions. The DPD predictions for the mass dependence of the diffusion coefficient for a single CG atactic polystyrene chain in toluene compare well with the experimental observations. The simulation results follow the power law dependence of the experimental diffusion coefficient on molar mass, which works well for medium to long chains. The only small deviation is observed for short chains and is attributed to the effect of chemically different end segments. Proper account of this effect requires a more sophisticated model with different CG units in the polymer chain. We did not include that in the present case study, polystyrenes used for practical applications are typically long-chain polymers. It might be needed, though, to model oligomers or polymers with different repeat units. Further extensions of this methodology to complex systems and solutions include, in particular, oligomeric polyelectrolytes in both aqueous and nonaqueous solvents and mixtures, as well as ionic liquids. Proper description beyond simplified models of complex solvation structures formed by association of polyelectrolytes and counterions in different solvents presents a significant challenge to molecular and CG simulation techniques. Our approach enables contraction of the degrees of freedom of both solvent and counterions by accounting for them in an effective pair potential between CG beads of polyelectrolyte. These effective pair potentials between neutralized CG beads account for dielectric screening by polar solvent and Debye screening by ionic species, and also include associative effects in the system. The remaining effort thus reduces to effective potential DPD simulation of a fluid of neutral CG polymeric chains at a lower packing fraction with a relatively short-range CG effective pair potential replacing the standard conservative force of DPD. We have successfully performed a preliminary test of this approach on multisolvent oligomeric polyelectrolyte gelator systems,57 which will be a subject of a subsequent publication.

A. BASICS OF DPD In DPD, to obtain the correct hydrodynamic behavior of the system, the resultant force Fi acting on a single bead is expressed as a sum of pairwise-additive conservative forces FCij , momentum-conserving dissipative forces FijD, and random forces FRij , where the sum runs over all neighboring beads within a certain cutoff distance Rc: Fi =

∑ (FCij + FijD + FijR )riĵ (A1)

j≠i

where r̂ij = rij/|rij| is a unit vector in the direction of rij. Hydrodynamics of the DPD fluid can be described by the Navier−Stokes equation, provided conservative forces are soft, and dissipative and random forces are short-ranged. For simple fluids without polymer chains, conservative forces FCij represent nonbonded interactions between DPD beads:

FCij = aijw(rij)riĵ

(A2)

specified with the weight function I


The Journal of Physical Chemistry B ⎧ ⎪1 − rij / R c for rij ≤ R c w(rij) = ⎨ ⎪ 0 for rij > R c ⎩

Article

kBTχ = 2α(aijuv − aijvv )(ρv + ρu )

where ρv and ρu are the number densities of solvent and solute beads. These densities can be used to define the cutoff distance Rc through the relation (ρv + ρu)R3c = N̅ , where the number N̅ is specified from practical reasons.57,70,71 The actual way of calculating the chemical potential difference μv(νu) − μv(0) is a matter of choice and availability. For example, in our previous studies on multisolvent oligomeric polyelectrolyte gelators,57,70,71 this was done by using the RISM-KH molecular theory of solvation.55 This approach is quite accurate and affordable in terms of computational resources. At the same time, one has to note that for most systems the requirement of equal volume of solute and solvent CG beads, which is a direct consequence of the lattice theory of polymer solutions,122 cannot be satisfied with mathematical exactness and is undesirable because it significantly restricts the freedom of coarse-graining. Overcoming this limitation has been our major motivation for the present study. For complex fluids with flexible polymer chains,75 besides nonbonded forces, the conservative force FC should also contain the contributions of bond stretching Fstr and bending Fbend:

(A3)

where aij > 0 is the repulsive force parameter between beads i and j at separation rij. The cutoff distance Rc can be obtained from the volume of the DPD bead as Rc = (ρ*VB/NA)1/3, where ρ* is the reduced density of DPD bead, VB is the molar volume of the bead, and NA is the Avogadro number. To find the parameters avv ij between solvent beads, the compressibility of the DPD solvent is typically matched to the value attributed to the real system, which is either calculated by some high-level theory or obtained experimentally.64,74 The approximate equation for pressure P of DPD fluid is given by P = ρkBT + αaρ2

(A4)

where kB is the Boltzmann constant, T is the thermodynamic temperature, α = 0.101 ± 0.001 is a characteristic parameter of DPD, a is a repulsive force parameter, and ρ is the bead number density per cubic cell Rc3. The isothermal dimensionless compressibility κT is given by64 κT−1 =

aρ ⎞ 1 ⎛ ∂P ⎞ ⎛ ∂ρ ⎞ 1 ⎛ ⎜1 + 2α ⎟ ⎜ ⎟ ⎜ ⎟ = kBT ⎝ ∂ρ ⎠T ⎝ ∂n ⎠T Nm ⎝ kBT ⎠

Fstr = −K str(rij − S0)riĵ

(A5)

where Nm is the number of molecules per bead, and n = ρNm is the number density of actual molecules per cubic cell R3c . The two equations above provide the necessary relationship between the mesoscopic model parameters and the compressibility of the system. Therefore, with κT determined from highlevel theory or experiment, a proper value of the conservative force coefficient avv ij can be obtained. With the requirement that beads of all of the components (both solute and solvent) have equal volume, this repulsion parameter should be set the same uu for all liquid components, that is, avv ij = aij  a, while the uv repulsion parameter aij between solute and solvent beads is determined from the Flory−Huggins theory.122 In particular, the excess free energy of mixing ΔAm is expressed in terms of the so-called mutual solubility parameter χ as ΔA m = Nv ln νv + Nu ln νu + χM 0νvνu kBT

(A9)

Fbend = −∇[Kbend(θ − θ0)2 ]

(A10)

where Kstr is the elastic constant of bond stretching and S 0 is the bond length at the bonding energy minimum; Kbend is the elastic constant of bending, θ is the angle between two adjacent bonds rij and rjk defined as cos θ = rij·rjk/rijrjk, and θ0 is the angle value at the bending potential minimum. These parameters are determined from higher-level (atomistic) theories. With the harmonic angle potential ϕθ in eq A10, the angle force is written as Fibend = −∇ϕθ = −

∂ϕθ ∂ cos θ ∂ cos θ ∂ri

Finally, the dissipative force defined as

(A6)

where Nv and Nu are the numbers of solvent and solute beads, respectively, νv and νu are the solvent and solute volume fractions with the normalization condition νv + νu = 1, and M0 is the total number of beads in the system. The parameter χ can be obtained experimentally, for example, as in ref 123. Yet more desirable for modeling is to get it from theory or simulation. With the free energy of mixing ΔAm available in terms of the chemical potential difference, the solubility parameter χ can be deduced from the following equation:122

FDij

(A11)

and the random force FRij are

FijD = −γw 2(rij)(vij·riĵ )riĵ

(A12)

FijR = σw(rij)ξijΔt −1/2 riĵ

(A13)

FDij

The dissipative force depends on both relative positions rij = rj−ri and relative velocities vij = vj − vi of the beads. The coefficient γ controlling the dissipative force strength is usually set equal to 1−10. The larger value produces a system that is very responsive to fluctuations in temperature and relaxes very quickly. The noise amplitude parameter σ determines the magnitude of the random pair force between the beads, ξij is the Gaussian distributed random variable with zero mean and unit variance, and Δt is the integration time step. Dissipative force FijD and random force FijR act as heat sink and source, respectively, so that their combined effect is a DPD thermostat.63,64 To make this particular thermostat conserve the momentum, the following relationship has to be obeyed:

μ (νu) − μv (0) 1 ⎛ ∂ΔA m ⎞ ⎜ ⎟ = v kBT ⎝ ∂Nv ⎠ kBT V

⎛ 1⎞ = ln νv + ⎜1 − ⎟ ln νu + χνu2 ⎝ M⎠

(A8)

(A7)

where μv(0) and μv(νu) are the chemical potentials of solvent molecules in pure solvent and in solution of polymer at a finite concentration, and M is the number of repeat units in each polymer chain. With the parameter χ available from eq A7, the interaction parameter auv ij is then obtained from the equation:

σ 2 = 2kBTγ

(A14)

As in any integration scheme, of great practical importance is a choice of Δt. A large value may cause beads to frequently change their nearest-neighbor assignments and destabilize the J



Article

system dynamics. On the other hand, a small value increases integration accuracy but requires much more time to perform the simulation. There have been studies into the effect of using various values for Δt; however, adaptive algorithms using values between 0.005 and 0.01 (in DPD time units) are sufficient for most practical purposes. Details on time unit conversion can be seen in our recent study.124

ϕtors(φ) = −kBT ln Ptors(φ)

which together constitute the bonded CG interaction potential: eff ϕbond (S , θ , φ) = ϕstr (S)ϕbend(θ )ϕtors(φ)

hαγ(r ) − ζαγ(r ) = [ωαμ(r ) + ζαμ(r )ρμ ]*cμν(r ) *[ωνγ(r ) + ρν hνγ(r )]

∫0

∞

gαγ(r ) = 1 + hαγ(r )

(B2a)

sin θ dθ Pbend(θ ) = 1

(B2b)

cαγ(r ) ≈ −uαγ(r )/(kBT )

ωαγ(k) = j0 (k Sαγ) for k Sαγ ≠ 0 = 0 for α ∈ a and γ ∈ c on a ≠ c

π

Each of the intrachain probability distribution functions sampled in the MD simulation and modeled with best multiGaussian fit corresponds through the Boltzmann inversion to a generalized intrachain PMF: (B3a)

ϕbend(θ ) = −kBT ln Pbend(θ )

(B3b)

(C4)

where j0(x) is the zeroth-order spherical Bessel function. Note that the DRISM integral eq C1 can be extended to flexible solute monomers and solvent molecules similarly to polymer RISM (PRISM) theory130−132 by parametrizing the corresponding intramolecular correlation functions ωαγ(r) in a multiGaussian form129 fitted to data from quantum chemical calculations, molecular simulations with classical force fields, or experiment. Both the total correlation function hαγ(r) and the intramolecular correlation function ω αγ (r) in eq C1 are

(B2c)

ϕstr(S) = −kBT ln Pstr(S)

(C3)

ωαγ(r) is the intramolecular matrix representing the geometry of molecular species, or the intramolecular correlation function between sites α, γ ∈ a of same molecular species of sort a normalized with ∫ dr ωαγ(r) = 1, and ωαγ = 0 between α ∈ a and γ ∈ c of different molecular species a ≠ c; ρα is the number density of interaction site α (of species a); and the asterisk “*” denotes convolution in the direct r-space and summation over repeating site indices. As mentioned in Appendix B, intramolecular degrees of freedom are treated as effectively separate from intermolecular ones. With both monomers and solvent molecules being relatively stiff, their geometry is adequately represented with rigid site separations S αγ fixed to the bond lengths at the energy minimum. The intramolecular matrix is then expressed in the direct space in terms of δ-functions, ωαγ(r) = δ(r − S αγ)/ (4πS 2αγ), and is specified in reciprocal k-space for numerical calculations as

π

∫−π dφ Ptors(φ) = 1

(C2)

measured in scattering experiments or sampled in molecular simulations; cαγ(r) is the site−site direct correlation function, which has the asymptotics of the site−site nonbonding interaction potential uαγ(r) (molecular force field):

(B1)

dS Pstr(S) = 1

(C1)

where hαγ(r) is the total correlation function at separation r between interaction sites α and γ enumerating all interaction sites on all sorts of molecular species in the solution, which is related to the site−site RDF:

To smooth MD simulation noise, the intrachain probability distribution functions Pstr(S ), Pbend(θ), and Ptors(φ) are modeled using a multi-Gaussian decomposition with expansion coefficients determined by nonlinear minimization of the average deviations.129 Each of the probability distribution functions is then normalized over its domain to unity:

∫0

(B4)

C. DRISM-KH MOLECULAR THEORY OF SOLVATION The site−site RDFs for the all-atom model of a solution system, which in general comprises solute monomers, molecular solvent or solvent mixture, counterions, and electrolyte, are obtained by converging the integral equations of the dielectrically consistent reference interaction site model (DRISM)52,53 complemented with the Kovalenko−Hirata (KH) closure relation.54−56 The DRISM integral equation52,53 can be cast in the form:

B. TREATMENT OF BONDED INTERACTIONS IN A CG POLYMER CHAIN There have been several previous studies114,125−127 of the mapping between the microscopic and mesoscopic models of polystyrene by using a bonded CG interaction potential ϕeff bond specifically constructed to address the inherent stiffness of degrees of freedom in CG structures.128 The net effect of all interactions in the system except for the bonding between CG polymer beads is represented with the nonbonded bead−bead CG effective pair potential ϕeff nb defined in section 2.2, which depends on the thermodynamic conditions and composition (temperature, pressure, concentration, pH, etc.). On the basis of the notion that bonded degrees of freedom in a CG chain are sufficiently stiff to prevent close contact between next and next−next neighbors, the interaction between beads within any four-unit fragment of the CG chain is modeled with a bonded CG interaction potential ϕeff bond sampled in a relatively short standard MD simulation of the all-atom system and used afterward. Solvent-mediated and multibody effects represented with nonbonded ϕeff nb are considered weak as compared to bonded interactions and thus are taken into account for bonded beads sufficiently well in the sampling of ϕeff bond. Furthermore, motion modes associated with intrachain degrees of freedom of bond stretching, bending, and torsion are treated not only as effectively separate from nonbonded interactions but also as mutually independent.128 Therefore, the probability distributions Pstr(S ) of bond length S , Pbend(θ) of bond angle θ, and Ptors(φ) of dihedral angle φ measured between interaction sites used for coarse-graining are treated as independent and uncorrelated: P(S , θ , φ) = Pstr(S)Pbend(θ )Ptors(φ)

(B3c)

K



Article

⎡ uαγ(r ) ⎤ + hαγ(r ) − cαγ(r )⎥ for gαγ(r ) ≤ 1 gαγ(r ) = exp⎢ − ⎣ kBT ⎦

renormalized with a dielectric bridge correction in an analytical form that ensures the given value and consistency of the dielectric constant obtained for solvent−solvent, solvent−ion, and ion−ion effective interactions in electrolyte solution.52,53 The renormalized dielectric correction is specified in the reciprocal space as52,53

=1−

uαγ(r ) kBT

+ hαγ(r ) − cαγ(r )

for gαγ(r ) > 1 (C10)

ζαγ(k) = j0 (kxα)j0 (kyα )j1 (kzα)hcj0 (kxγ )j0 (kyγ )j1 (kz γ ), for α ,

where the site−site interaction potential uαγ(r) (all-atom force field) is modeled typically with Coulomb and Lennard-Jones terms. The KH closure C10 couples in a nontrivial way the socalled hypernetted chain (HNC) approximation and the meanspherical approximation (MSA) linearization. The former is held for spatial regions of density depletion gαγ(r) < 1 (including the repulsive core), and the latter is automatically applied to spatial regions of solvent density enrichment gαγ(r) > 1 such as strong peaks of association and long-range tails of near-critical fluid phases, while keeping the long-range asymptotics peculiar to both the HNC and the MSA. (The distribution function and its first derivative are continuous at the joint boundary gαγ(r) = 1 by construct.) The DRISM-KH eqs C1 and C10 thus possess the same dielectrically consistent asymptotics (eq C3) as the originally derived DRISM-HNC theory,52,53 but extend the description to solution systems with complex associative species in a wide range of composition and thermodynamic conditions not amenable to HNC. The latter overestimates associative attraction and nanostructure density inhomogeneities and usually diverges for such systems. As distinct, the KH closure consistently accounts for both electrostatic and nonpolar features of solvation (associative and steric effects) in simple and complex liquids, ionic liquids, and nonelectrolyte and electrolyte solutions in various chemical systems.56 In particular, the KH closure approximation successfully reproduced a number of experimentally observed phenomena in soft matter systems and liquid interfaces, including the structural transitions and related thermodynamic anomalies for the formation of micromicelles and the transition from the tetrahedral to zigzag hydrogen-bonding network in water−alcohol mixtures in the whole range of concentrations,60,61 the microscopic structure of gels formed by oligomeric polyelectrolyte gelators in different solvents,57 phase equilibria in bulk and nanoconfined fluids,134,135 and microscopic structure of interfaces of nonpolar and polar molecular liquids.136,137 The KH closure underestimates the height of strong associative peaks of the distribution functions because of the MSA linearization applied to them.54,138 However, it somewhat widens the peaks and so quite accurately reproduces the coordination numbers of solvation structure. This has been demonstrated for different systems, including microheterogeneities (micromicellar aggregates) in water−alcohol solutions, 56,60 and in the context of the 3D-RISM-KH theory54−56,94 for solvation shells of a metal−water54,55 and metal oxide−water139,140 interfaces, and structural water solvent localized in biomolecular confinement.58,59 For example, the coordination numbers of water strongly bound to the MgO surface are calculated from the 3D-RISM-KH theory with a 90% accuracy and the peak positions within a 0.5 Å deviation, as compared to MD simulations.139 The 3D solvation map of function-related structural water in the GroEL chaperon complex (shown to be strongly correlated to the rate of protein folding inside the chaperon cavity) obtained in an expensive MD simulation with explicit solvent involving ∼1

γ ∈ polar a = 0 otherwise (C5)

where j0(x) and j1(x) are, respectively, the zeroth- and firstorder spherical Bessel functions, and rα = (xα, yα, zα) and rγ = (xγ, yγ, zγ) are the Cartesian coordinates of partial site charges qα and qγ on the same polar species a with respect to its molecular origin, such that its dipole moment da∑α∈aqαrα is oriented along the z-axis, da = (0, 0, da). Note that the renormalized dielectric correction (eq C5) is nonzero only for solution species that possess a dipole moment and thus are responsible for the dielectric response in the DRISM approach, and ζαγ = 0 for α ∈ a and γ ∈ c on nonpolar or different species a ≠ c. The value of the envelope function hc(k) at k = 0 determines the dielectric constant of the solution, and the rest of the function is assumed in a smooth nonoscillatory form quickly falling off at wavevectors k larger than an inverse of characteristic size l of liquid molecules: hc(k) = A exp[−(lk)2 /4]

(C6)

where the amplitude A is related to the solution dielectric constant ε, which is a phenomenological parameter specified at input: A=

⎞ 1 ⎛ε ⎜ − 3⎟ ρpolar ⎝ y ⎠

(C7)

The form C5 can be extended to mixed solvents133 by using the total number density of solution polar species:

ρpolar =

∑

ρa (C8)

a ∈ polar

and the mixture dielectric susceptibility: y=

4π 9kBT

∑ a ∈ polar

ρa da2 (C9)

The parameter l in the envelope function C6 specifies the characteristic separation from a molecule in solution below which the dielectric correction C5 is switched off so as not to distort the short-range solvation structure. It can be set to about l = 1 Å for water solvent; however, in solvent of larger molecules such as for example octanol or in the presence of such cosolvent, it should be increased to about l = 10 Å so as to avoid “ghost” associative peaks appearing in the radial distributions if the dielectric correction C5 interferes with the intramolecular structure of the large solvent species. The DRISM integral eq C1 is complemented with the KH closure,54−56 which in the site−site radial version reads L



■

million atoms is reproduced from the 3D-RISM-KH theory in a relatively short calculation on a workstation with an accuracy of over 90% correlation for the 3D density map and about 98% correlation for the 3D positions of density maxima.59 In numerical implementation, special treatment is required for electrostatic forces in solution with polar molecular solvent and ionic species. We separate out and handle analytically the electrostatic asymptotics in the direct and reciprocal space from all of the correlation functions in the DRISM-KH eqs C1 and C10.55,56 The remaining short-range terms are discretized on a uniform radial grid with resolution 0.01−0.1 Å, and their convolution is handled with the technique of sine fast Fourier transform, or discrete sine transform. The DRISM-KH integral eqs C1 and C10 are solved for the site−site direct correlation functions cαγ(r) of the all-atom solution of monomers by using the modified direct inversion in iterated space (MDIIS) accelerated numerical solver,55,93,94 which also yields the site−site RDFs gαγ(r). The thermodynamic properties of the system are then analytically obtained in terms of a simple integral of the correlation functions.55,56

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Attinger, S., Koumoutsakos, P., Eds. Multiscale Modelling and Simulation; Lecture Notes in Computational Science and Engineering; Springer: Berlin, 2004; Vol. 39. (2) Liu, W. K.; Karpov, E. G.; Park, H. S. Nano Mechanics and Materials: Theory, Multiscale Methods and Applications; John Wiley & Sons, Ltd.: Chichester, 2006. (3) Guo, Z. X., Ed. Multiscale Materials Modelling. Fundamentals and Applications; Woodhead Publishing Ltd.: Cambridge, 2007. (4) Kwon, Y. W., Allen, D. H., Talreja, R., Eds. Multiscale Modeling and Simulation of Composite Materials and Structures; Springer: New York, 2008. (5) Ross, R. B., Mohanty, S., Eds. Multiscale Simulation Methods for Nanomaterials; John Wiley & Sons, Inc.: Hoboken, 2008. (6) Steinhauser, M. O. Computational Multiscale Modeling of Fluids and Solids: Theory and Applications; Springer: Berlin, 2008. (7) Voth, G. A., Ed. Coarse-Graining of Condensed Phase and Biomolecular Systems; CRC Press: Boca Raton, FL, 2008. (8) Engquist, B., Lötstedt, P., Runborg, O., Eds. Multiscale Modeling and Simulation in Science; Lecture Notes in Computational Science and Engineering; Springer: Berlin, 2009; Vol. 66. (9) Fish, J., Ed. Multiscale Methods Bridging the Scales in Science and Engineering; Oxford University Press: Oxford, 2009. (10) King, M. R., Gee, D. J., Eds. Multiscale Modeling of Particle Interactions. Applications in Biology and Nanotechnology; John Wiley & Sons, Inc.: Hoboken, 2010. (11) Derosa, P., Cagin, T., Eds. Multiscale Modeling: From Atoms to Devices; CRC Press: Boca Raton, FL, 2011. (12) E, W. Principles of Multiscale Modeling; Cambridge University Press: Cambridge, 2011. (13) Li, J.; Ge, W.; Wang, W.; Yang, N.; Liu, X.; Wang, L.; He, X.; Wang, X.; Wang, J.; Kwauk, M. From Multiscale Modeling to MesoScience: A Chemical Engineering Perspective. Principles, Modeling, Simulation, and Application; Springer: Berlin, 2013. (14) Li, S., Qian, D., Eds. Multiscale Simulations and Mechanics of Biological Materials; John Wiley & Sons Ltd.: Chichester, 2013. (15) Monticelli, L., Salonen, E., Eds. Biomolecular Simulations: Methods and Protocols. Methods in Molecular Biology; Humana Press: New York, 2013; Vol. 924. (16) Müller-Plathe, F. Coarse-Graining in Polymer Simulation: From the Atomistic to the Mesoscopic Scale and Back. ChemPhysChem 2002, 3, 754−769. (17) Karakasidis, T. E.; Charitidis, C. A. Multiscale Modeling in Nanomaterials Science. Mater. Sci. Eng., C 2007, 27, 1082−1089. (18) Ortoleva, P.; Singharoy, A.; Pankavich, S. Hierarchical Multiscale Modeling of Macromolecules and Their Assemblies. Soft Matter 2013, 9, 4319−4335. (19) Salmon, P. S.; Zeidler, A. Identifying and Characterising the Different Structural Length Scales in Liquids and Glasses: An Experimental Approach. Phys. Chem. Chem. Phys. 2013, 15, 15286− 15308. (20) Yip, S.; Short, M. P. Multiscale Materials Modelling at the Mesoscale. Nat. Mater. 2013, 12, 774−777. (21) Glinski, G.; Bailey, C.; Pericleous, K. A. Application of Coupled Continuum-Mesoscopic Computational Methods for the Simulation of Complex Fluids in Industrial Processes. Computational Modeling of Materials, Minerals and Metals Processing; Conference on Computational Modeling of Materials, Minerals and Metals Processing: San Diego, CA, Sep. 23−26, 2001; pp 149−160. (22) McGrother, S.; Ming, L. Y.; Goldbeck-Wood, G. Mesoscale Simulations: Industrial Applications. In Mesoscale Phenomena in Fluid Systems; Case, F., Alexandridis, P., Eds.; American Chemical Society: Washington, DC, 2003; Chapter 15, pp 227−241. (23) Fermeglia, M.; Pricl, S. Multiscale Modeling for Polymer Systems of Industrial Interest. Prog. Org. Coat. 2007, 58, 187−199. (24) Fermeglia, M.; Pricl, S. Multiscale Molecular Modeling in Nanostructured Material Design and Process System Engineering. Comput. Chem. Eng. 2009, 33, 1701−1710.

D. FLOWCHART DIAGRAM FOR DPD WITH CG EFFECTIVE PAIR POTENTIAL FROM DRISM-KH

■

Article

ACKNOWLEDGMENTS

This work was supported by the National Institute for Nanotechnology and the University of Alberta. The high performance computing resources for this work were provided by WestGrid − Compute/Calcul Canada. Figure 1 was created by using Materials Studio Visualizer,31 and Figures 3−6 were produced using Origin.121 We express our gratitude to Dr. Xiangjun Liu for his commitment to the project in its early stage. M



Article

(25) Gillet, G.; Vitrac, O.; Desobry, S. Prediction of Partition Coefficients of Plastic Additives Between Packaging Materials and Food Simulants. Ind. Eng. Chem. Res. 2010, 49, 7263−7280. (26) Spaeth, J. R.; Dale, T.; Kevrekidis, I. G.; Panagiotopoulos, A. Z. Coarse-Graining of Chain Models in Dissipative Particle Dynamics Simulations. Ind. Eng. Chem. Res. 2011, 50, 69−77. (27) Markutsya, S.; Fox, R. O.; Subramaniam, S. Coarse-Graining Approach to Infer Mesoscale Interaction Potentials from Atomistic Interactions for Aggregating Systems. Ind. Eng. Chem. Res. 2012, 51, 16116−16134. (28) Morohoshi, K. Innovative High-Tech Process of Polymer Design: Application to Fuel Cell Polyelectrolyte Membrane. Kobunshi Ronbunshu 2012, 69, 493−501 (in Japanese). (29) Chen, C.; Lu, K.; Li, X.; Dong, J.; Lu, J.; Zhuang, L. A ManyBody Dissipative Particle Dynamics Study of Fluid-Fluid Spontaneous Capillary Displacement. RSC Adv. 2014, 4, 6545−6555. (30) Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G. L.; Cococcioni, M.; Dabo, I.; et al. QUANTUM ESPRESSO: A Modular and Open-Source Software Project for Quantum Simulations of Materials. J. Phys.: Condens. Matter 2009, 21, 395502 http://www.quantum-espresso.org. (31) Accelrys Software Inc., Materials Studio Modeling Environment, Release 5.5; Accelrys Software Inc: San Diego, CA, 2010; http://www. accelrys.com. (32) Binder, K.; Paul, W.; Santos, S.; Suter, U. W. Coarse-Graining Techniques. In Simulation Methods for Polymers; Kotelyanskii, M., Theodorou, D. N., Eds.; Marcel Dekker: New York, 2004; Chapter 15, pp 491−510. (33) Peter, C.; Kremer, K. Multiscale Simulation of Soft Matter Systems − from the Atomistic to the Coarse-Grained Level and Back. Soft Matter 2009, 5, 4357−4366. (34) Morrow, T. Overcoming Large Time- and Length-Scale Challenges in Molecular Modeling. In Multiscale Modeling: From Atoms to Devices; Derosa, P., Cagin, T., Eds.; CRC Press: Boca Raton, FL, 2011; Chapter 1, pp 1−12. (35) Karimi-Varzaneh, H.; Mü ller-Plathe, F. Coarse-Grained Modeling for Macromolecular Chemistry. In Multiscale Molecular Methods in Applied Chemistry; Kirchner, B., Vrabec, J., Eds.; Topics in Current Chemistry; Springer: Berlin, 2012; Vol. 307, pp 295−321. (36) Mazenko, G. F. Coarse Graining and Effective Hamiltonians. Fluctuations, Order, and Defects; John Wiley & Sons: Hoboken, 2003; Chapter 2, pp 57−105. (37) Mazenko, G. F. Coarse Graining and Effective Hamiltonians, and the Renormalization Group. Fluctuations, Order, and Defects; John Wiley & Sons: Hoboken, 2003; Chapter 3, pp 106−154. (38) Louis, A. A. Beware of Density Dependent Pair Potentials. J. Phys.: Condens. Matter 2002, 14, 9187. (39) D’Adamo, G.; Pelissetto, A.; Pierleoni, C. Coarse-Graining Strategies in Polymer Solutions. Soft Matter 2012, 8, 5151−5167. (40) D’Adamo, G.; Pelissetto, A.; Pierleoni, C. Predicting the Thermodynamics by Using State-Dependent Interactions. J. Chem. Phys. 2013, 138, 234107. (41) Ercolessi, F.; Adams, J. B. Interatomic Potentials from FirstPrinciples Calculations: The Force-Matching Method. Europhys. Lett. 1994, 26, 583−588. (42) Izvekov, S.; Voth, G. A. A Multiscale Coarse-Graining Method for Biomolecular Systems. J. Phys. Chem. B 2005, 109, 2469−2473. (43) Izvekov, S.; Voth, G. A. Multiscale Coarse Graining of LiquidState Systems. J. Chem. Phys. 2005, 123, 134105. (44) Lu, L.; Dama, J. F.; Voth, G. A. Fitting Coarse-Grained Distribution Functions through an Iterative Force-Matching Method. J. Chem. Phys. 2013, 139, 121906. (45) Snook, I. The Langevin and Generalized Langevin Approach to the Dynamics of Atomic, Polymeric and Colloidal Systems; Elsevier: Amsterdam, 2007. (46) Boon, J. P.; Yip, S. Molecular Hydrodynamics; Wiley Interscience: New York, 1980; [Reprinted by Dover Publications, New York, 1991]. (47) Hansen, J. P.; McDonald, I. Theory of Simple Liquids, 3rd ed.; Elsevier: Amsterdam, 2006.

(48) Bolhuis, P. G.; Louis, A. A. How to Derive and Parameterize Effective Potentials in Colloid-Polymer Mixtures. Macromolecules 2002, 35, 1860−1869. (49) Bolhuis, P. G.; Louis, A. A.; Hansen, J.-P. Many-Body Interactions and Correlations in Coarse-Grained Descriptions of Polymer Solutions. Phys. Rev. E 2001, 64, 021801. (50) Hirata, F. Theory of Molecular Liquids. In Molecular Theory of Solvation; Hirata, F., Ed.; Understanding Chemical Reactivity; Kluwer: Dordrecht, 2003; Vol. 24, Chapter 1, pp 1−60. (51) Kast, S. M.; Hauptmann, W.; Schilling, B. Design of Reduced Molecular Models by Integral Equation Theory. Chem. Phys. Lett. 2006, 418, 202−207. (52) Perkyns, J. S.; Pettitt, B. M. A Dielectrically Consistent Interaction Site Theory for Solvent-Electrolyte Mixtures. Chem. Phys. Lett. 1992, 190, 626−630. (53) Perkyns, J. S.; Pettitt, B. M. A Site-Site Theory for Finite Concentration Saline Solutions. J. Chem. Phys. 1992, 97, 7656−7666. (54) Kovalenko, A.; Hirata, F. Self-Consistent Description of a MetalWater Interface by the Kohn-Sham Density Functional Theory and the Three-Dimensional Reference Interaction Site Model. J. Chem. Phys. 1999, 110, 10095−10112. (55) Kovalenko, A. Three-Dimensional RISM Theory for Molecular Liquids and Solid-Liquid Interfaces. In Molecular Theory of Solvation; Hirata, F., Ed.; Understanding Chemical Reactivity; Kluwer: Dordrecht, 2003; Vol. 24; Chapter 4, pp 169−275. (56) Kovalenko, A. Multiscale Modeling of Solvation in Chemical and Biological Nanosystems and in Nanoporous Materials. Pure Appl. Chem. 2013, 85, 159−199. (57) Kovalenko, A.; Kobryn, A. E.; Gusarov, S.; Lyubimova, O.; Liu, X.; Blinov, N.; Yoshida, M. Molecular Theory of Solvation for Supramolecules and Soft Matter Structures: Application to Ligand Binding, Ion Channels, and Oligomeric Polyelectrolyte Gelators. Soft Matter 2012, 8, 1508−1520. (58) Yoshida, N.; Imai, T.; Phongphanphanee, S.; Kovalenko, A.; Hirata, F. Molecular Recognition in Biomolecules Studied by Statistical-Mechanical Integral-Equation Theory of Liquids. J. Phys. Chem. B 2009, 113, 873−886. (59) Stumpe, M. C.; Blinov, N.; Wishart, D.; Kovalenko, A.; Pande, V. S. Calculation of Local Water Densities in Biological Systems: A Comparison of Molecular Dynamics Simulations and the 3D-RISMKH Molecular Theory of Solvation. J. Phys. Chem. B 2011, 115, 319− 328. (60) Yoshida, K.; Yamaguchi, T.; Kovalenko, A.; Hirata, F. Structure of Tert-Butyl Alcohol-Water Mixtures Studied by the RISM Theory. J. Phys. Chem. B 2002, 106, 5042−5049. (61) Omelyan, I.; Kovalenko, A.; Hirata, F. Compressibility of TertButyl Alcohol-Water Mixtures: The RISM Theory. J. Theor. Comput. Chem. 2003, 2, 193−203. (62) Hoogerbrugge, P. J.; Koelman, J. M. V. A. Simulating Microscopic Hydrodynamic Phenomena with Dissipative Particle Dynamics. Europhys. Lett. 1992, 19, 155−160. (63) Espanol, P.; Warren, P. Statistical Mechanics of Dissipative Particle Dynamics. Europhys. Lett. 1995, 30, 191−196. (64) Groot, R. D.; Warren, P. B. Dissipative Particle Dynamics: Bridging the Gap Between Atomistic and Mesoscopic Simulation. J. Chem. Phys. 1997, 107, 4423−4435. (65) Sims, J. S.; Martys, N. S. Simulation of Sheared Suspensions with a Parallel Implementation of QDPD. J. Res. Natl. Inst. Stand. Technol. 2004, 109, 267−277. (66) Martys, N. S. Study of a Dissipative Particle Dynamics Based Approach for Modeling Suspensions. J. Rheol. 2005, 49, 401−424. (67) Goujon, F.; Malfreyt, P. J.; Tildesley, D. Mesoscopic Simulation of Entangled Polymer Brushes under Shear: Compression and Rheological Properties. Macromolecules 2009, 42, 4310−4318. (68) Levine, Y. K.; Gomes, A. E.; Martines, A. F.; Polimeno, A. A Dissipative Particle Dynamics Description of Liquid-Crystalline Phases. I. Methodology and Applications. J. Chem. Phys. 2005, 122, 144902. N



Article

molecular Chemistry; Gokel, G. W., Ed.; JAI Press Inc.: Stamford, CT, 2000; Vol. 6, Chapter 2, pp 41−156. (89) Qin, J.; Milner, S. T.; Stephanou, P. S.; Mavrantzas, V. G. Effects of Tube Persistence Length on Dynamics of Mildly Entangled Polymers. J. Rheol. 2012, 56, 707−723. (90) Lyubartsev, A. P.; Laaksonen, A. Calculation of Effective Interaction Potentials from Radial Distribution Functions: A Reverse Monte Carlo Approach. Phys. Rev. E 1995, 52, 3730−3737. (91) Lyubartsev, A. P.; Laaksonen, A. Effective Potentials for IonDNA Interactions. J. Chem. Phys. 1999, 111, 11207−11215. (92) Mirzoev, A.; Lyubartsev, A. P. Effective Solvent Mediated Potentials of Na+ and Cl− Ions in Aqueous Solution: Temperature Dependence. Phys. Chem. Chem. Phys. 2011, 13, 5722−5727. (93) Kovalenko, A.; Ten-no, S.; Hirata, F. Solution of the ThreeDimensional RISM/HNC Equations for SPC Water by the Modified Method of Direct Inversion in the Iterative Subspace. J. Comput. Chem. 1999, 20, 928−936. (94) Kovalenko, A.; Hirata, F. Potentials of Mean Force of Simple Ions in Ambient Aqueous Solution. I. Three-Dimensional Reference Interaction Site Model Approach. J. Chem. Phys. 2000, 112, 10391− 10402. (95) Urakawa, O. Self-Diffusion and Viscosity of Low Molecular Weight Polystyrene Over a Wide Temperature Range. Macromolecules 2004, 37, 1558−1564. (96) Harmandaris, V. A.; Kremer, K. Dynamics of Polystyrene Melts through Hierarchical Multiscale Simulation. Macromolecules 2009, 42, 791−802. (97) Harmandaris, V. A.; Floudas, G.; Kremer, K. Temperature and Pressure Dependence of Polystyrene Dynamics through Molecular Dynamics Simulations and Experiments. Macromolecules 2011, 44, 393−402. (98) Karimi-Varzaneh, H. A.; van der Vegt, N. F. A.; Müller-Plathe, F.; Carbone, P. How Good are Coarse-Grained Polymer Models? A Comparison for Atactic Polystyrene. ChemPhysChem 2012, 13, 3428− 3439. (99) Sun, Q.; Faller, R. Systematic Coarse-Graining of Atomistic Models for Simulation of Polymeric Systems. Comput. Chem. Eng. 2005, 29, 2380−2385. (100) Rieger, J. The Glass Transition Temperature of Polystyrene. J. Therm. Anal. 1996, 46, 965−972. (101) Pronk, S.; Páll, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D.; et al. GROMACS 4.5: A High-Throughput and Highly Parallel Open Source Molecular Simulation Toolkit. Bioinformatics 2013, 29, 845−854. (102) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435− 447. (103) van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, Flexible, and Free. J. Comput. Chem. 2005, 26, 1701−1718. (104) Lindahl, E.; Hess, B.; van der Spoel, D. GROMACS 3.0: A Package for Molecular Simulation and Trajectory Analysis. J. Mol. Model. 2001, 7, 306−317. (105) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. GROMACS: A Message-Passing Parallel Molecular Dynamics Implementation. Comput. Phys. Commun. 1995, 91, 43−56. (106) Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling Through Velocity Rescaling. J. Chem. Phys. 2007, 126, 014101. (107) Nose, S.; Klein, M. L. Constant Pressure Molecular Dynamics for Molecular Systems. Mol. Phys. 1983, 50, 1055−1076. (108) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182−7190. (109) Jorgensen, W. L.; Madura, J. D.; Swenson, C. J. Optimized Intermolecular Potential Functions for Liquid Hydrocarbons. J. Am. Chem. Soc. 1984, 106, 6638−6646.

(69) Ibergay, C.; Malfreyt, P.; Tildesley, D. J. Mesoscale Modeling of Polyelectrolyte Brushes with Salt. J. Phys. Chem. B 2010, 114, 7274− 7284. (70) Liu, X.; Lyubimova, O.; Kobryn, A. E.; Gusarov, S.; Kovalenko, A. Mesoscopic Study of Dynamics and Gelation Ability of Oligomeric Electrolyte Gelator with Dissipative Particle Dynamics. Procedia Comput. Sci. 2011, 4, 1031−1038. (71) Lyubimova, O.; Liu, X.; Gusarov, S.; Kobryn, A. E.; Kovalenko, A. Solvation Structure and Gelation Ability of Polyelectrolytes: Predictions by Quantum Chemistry Methods and Integral Equation Theory of Molecular Liquids. Procedia Comput. Sci. 2011, 4, 1186− 1192. (72) Vishnyakov, A.; Talaga, D. S.; Neimark, A. V. DPD Simulation of Protein Conformations: From α-Helices to β-Structures. J. Phys. Chem. Lett. 2012, 3, 3081−3087. (73) Vishnyakov, A.; Lee, M.-T.; Neimark, A. V. Prediction of the Critical Micelle Concentration of Nonionic Surfactants by Dissipative Particle Dynamics Simulations. J. Phys. Chem. Lett. 2013, 4, 797−802. (74) Groot, R. D.; Rabone, K. L. Mesoscopic Simulation of Cell Membrane Damage, Morphology Change and Rupture by Nonionic Surfactants. Biophys. J. 2001, 81, 725−736. (75) de Meyer, F. J.; Venturoli, M.; Smit, B. Molecular Simulations of Lipid-Mediated Protein-Protein Interactions. Biophys. J. 2008, 95, 185−865. (76) Guigas, G.; Worozova, D.; Weiss, M. Exploring Membrane and Protein Dynamics with Dissipative Particle Dynamics. Adv. Protein Chem. Struct. Biol. 2011, 85, 143−182. (77) Yamamoto, S.; Maruyama, Y.; Hyodo, S. Dissipative Particle Dynamics Study of Spontaneous Vesicle Formation of Amphiphilic Molecules. J. Chem. Phys. 2002, 116, 5842−5849. (78) Ortiz, V.; Nielsen, S. O.; Discher, D. E.; Klein, M. L. Dissipative Particle Dynamics Simulations of Polymersomes. J. Phys. Chem. B 2005, 109, 17708−17714. (79) Dutt, M.; Kuksenok, O.; Nayhouse, M. J.; Little, S. R.; Balazs, A. C. Modeling the Self-Assembly of Lipids and Nanotubes in Solution: Forming Vesicles and Bicelles with Transmembrane Nanotube Channels. ACS Nano 2011, 5, 4769−4782. (80) Groot, R. D. Electrostatic Interactions in Dissipative Particle Dynamics-Simulation of Polyelectrolytes and Anionic Surfactants. J. Chem. Phys. 2003, 118, 11265−11277. (81) Chen, S.; Guo, C.; Hu, G.; Liu, H.; Liang, X.; Wang, J.; Ma, J.; Zheng, L. Dissipative Particle Dynamics Simulation of Gold Nanoparticles Stabilization by PEO-PPO-PEO Block Copolymer Micelles. Colloid Polym. Sci. 2007, 285, 1543−1552. (82) Zhong, C.; Liu, D. Understanding Multicompartment Micelles Using Dissipative Particle Dynamics Simulation. Macromol. Theor. Simul. 2007, 16, 141−157. (83) Otter, W. K. D.; Clarke, J. H. R. Simulation of Polymers by Dissipative Particle Dynamics. In Simulation Methods for Polymers; Kotelyanskii, M., Theodorou, D. N., Eds.; Marcel Dekker: New York, 2004; Chapter 17, pp 559−574. (84) Satoh, A. Practice of Dissipative Particle Dynamics Simulations. Introduction to Practice of Molecular Simulation: Molecular Dynamics, Monte Carlo, Brownian Dynamics, Lattice Boltzmann, Dissipative Particle Dynamics; Elsevier Insights; Elsevier: Tokyo, 2011; Chapter 6, pp 187−218. (85) Lu, Z.-Y.; Wang, Y.-L. An Introduction to Dissipative Particle Dynamics. In Biomolecular Simulations: Methods and Protocols; Monticelli, L., Salonen, E., Eds.; Methods in Molecular Biology; Humana Press: New York, 2013; Vol. 924, Chapter 24, pp 617−633. (86) Brini, E.; Algaer, E. A.; Ganguly, P.; Li, C.; Rodriguez-Ropero, F.; van der Vegt, N. F. A. Systematic Coarse-Graining Methods for Soft Matter Simulations − A Review. Soft Matter 2013, 9, 2108−2119. (87) Soddemann, T.; Dünweg, B.; Kremer, K. Dissipative Particle Dynamics: A Useful Thermostat for Equilibrium and Nonequilibrium Molecular Dynamics Simulations. Phys. Rev. E 2003, 68, 046702. (88) Feiters, M. C.; Nolte, R. J. M. Chiral Self-Assembled Structures from Biomolecules and Synthetic Analogues. In Advances in SupraO



Article

(110) Jorgensen, W. L.; Severance, D. L. Aromatic-Aromatic Interactions: Free Energy Profiles for the Benzene Dimer in Water, Chloroform, and Liquid Benzene. J. Am. Chem. Soc. 1990, 112, 4768− 4774. (111) Looyenga, H. Dielectric Constants of Heterogeneous Mixtures. Physica 1965, 31, 401−406. (112) LAMMPS. Large-scale Atomic/Molecular Massively Parallel Simulator, 2013; http://lammps.sandia.gov. (113) Humphrey, W.; Dalke, A.; Schulten, K. VMD − Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38, TopoTools by Axel Kohlmeyer. (114) Qian, H.-J.; Carbone, P.; Chen, X.; Karimi-Varzaneh, H. A.; Liew, C. C.; Müller Plathe, F. Temperature-Transferable CoarseGrained Potentials for Ethylbenzene, Polystyrene, and their Mixtures. Macromolecules 2008, 41, 9919−9929. (115) Farah, K.; Fogarty, A. C.; Böhm, M. C.; Müller-Plathe, F. Temperature Dependence of Coarse-Grained Potentials for Liquid Hexane. Phys. Chem. Chem. Phys. 2011, 13, 2894−2902. (116) Rauch, J.; Köhler, W. Collective and Thermal Diffusion in Dilute, Semidilute, and Concentrated Solutions of Polystyrene in Toluene. J. Chem. Phys. 2003, 119, 11977−11988. (117) Schäfer, R. Stochastische breitbandmodulierte transiente holographische Gitter zur Untersuchung des Ludwig-Soret-Effekts in mehrkomponentigen Flüssigkeiten. Ph.D. thesis, Johannes-Gutenberg Universität, Mainz, 1997. (118) Köhler, W.; Rosenauer, C.; Rossmanith, P. Holographic Grating Study of Mass and Thermal Diffusion of Polystyrene/Toluene Solutions. Int. J. Thermophys. 1995, 16, 11−21. (119) Rauch, J.; Köhler, W. On the Molar Mass Dependence of the Thermal Diffusion Coefficient of Polymer Solutions. Macromolecules 2005, 38, 3571−3573. (120) Bayramoglu, B.; Faller, R. Coarse-Grained Modeling of Polystyrene in Various Environments by Iterative Boltzmann Inversion. Macromolecules 2012, 45, 9205−9219. (121) Origin, Version 8.1.13.88, 2009; http://www.originlab.com. (122) Hill, T. L. An Introduction to Statistical Thermodynamics; Addison-Wesley Series in Chemistry; Addison-Wesley: Massachusetts, 1960; [Reprinted by Dover Publications, New York, 1986]. (123) Emerson, J. A.; Toolan, D. T. W.; Howse, J. R.; Furst, E. M.; Epps, T. H. Determination of Solvent-Polymer and Polymer-Polymer Flory-Huggins Interaction Parameters for Poly(3-Hexylthiophene) via Solvent Vapor Swelling. Macromolecules 2013, 46, 6533−6540. (124) Nikolić, D.; Moffat, K. A.; Farrugia, V. M.; Kobryn, A. E.; Gusarov, S.; Wosnick, J. H.; Kovalenko, A. Multi-Scale Modeling and Synthesis of Polyester Ionomers. Phys. Chem. Chem. Phys. 2013, 15, 6128−6138. (125) Harmandaris, V. A.; Reith, D.; van der Vegt, N. F. A.; Kremer, K. Comparison Between Coarse-Graining Models for Polymer Systems: Two Mapping Schemes for Polystyrene. Macromol. Chem. Phys. 2007, 208, 2109−2120. (126) Fritz, D.; Harmandaris, V. A.; Kremer, K.; van der Vegt, N. F. A. Coarse-Grained Polymer Melts Based on Isolated Atomistic Chains: Simulation of Polystyrene of Different Tacticities. Macromolecules 2009, 42, 7579−7588. (127) Rossi, G.; Monticelli, L.; Puisto, S. R.; Vattulainen, I.; AlaNissila, T. Coarse-Graining Polymers with the MARTINI Force-Field: Polystyrene as a Benchmark Case. Soft Matter 2011, 7, 698−708. (128) Li, C.; Shen, J.; Peter, C.; van der Vegt, N. F. A. A Chemically Accurate Implicit-Solvent Coarse-Grained Model for Polystyrenesulfonate Solutions. Macromolecules 2012, 45, 2551−2561. (129) Milano, G.; Goudeau, S.; Müller-Plathe, F. Multicentered Gaussian-Based Potentials for Coarse-Grained Polymer Simulations: Linking Atomistic and Mesoscopic Scales. J. Polym. Sci., Part B 2005, 43, 871−885. (130) Schweizer, K. S.; Curro, J. G. Integral-Equation Theory of the Structure of Polymer Melts. Phys. Rev. Lett. 1987, 58, 246−249. (131) Schweizer, K. S.; Curro, J. G. Integral Equation Theory of the Structure and Thermodynamics of Polymer Blends. J. Chem. Phys. 1989, 91, 5059.

(132) Schweizer, K. S.; Curro, J. G. Integral Equation Theories of the Structure, Thermodynamics, and Phase Transitions of Polymer Fluids. In Advances in Chemical Physics; Prigogine, I., Rice, S. A., Eds.; John Wiley & Sons, Inc.: New York, 1997; Vol. 98, pp 1−142. (133) Kvamme, B. Interaction-Site Representation of Polar Mixtures and Electrolyte Solutions. Int. J. Thermophys. 1995, 16, 743−750. (134) Kovalenko, A.; Hirata, F. First-Principles Realization of a van der Waals-Maxwell Theory for Water. Chem. Phys. Lett. 2001, 349, 496−502. (135) Kovalenko, A.; Hirata, F. Towards a Molecular Theory for the van der Waals-Maxwell Description of Fluid Phase Transitions. J. Theor. Comput. Chem. 2002, 1, 381−406. (136) Kovalenko, A.; Hirata, F. A Molecular Theory of Liquid Interfaces. Phys. Chem. Chem. Phys. 2005, 7, 1785−1793. (137) Kovalenko, A.; Hirata, F. A Molecular Theory of Solutions at Liquid Interfaces. In Interfacial Nanochemistry: Molecular Science and Engineering at Liquid-Liquid Interfaces; Watarai, H., Teramae, N., Sawada, T., Eds.; Nanostructure Science and Technology; Kluwer Academic/Plenum Publishers: New York, 2005; Chapter 5, pp 97− 125. (138) Perkyns, J. S.; Lynch, G. C.; Howard, J. J.; Pettitt, B. M. Protein Solvation from Theory and Simulation: Exact Treatment of Coulomb Interactions in Three-Dimensional Theories. J. Chem. Phys. 2010, 132, 064106. (139) Shapovalov, V.; Truong, T. N.; Kovalenko, A.; Hirata, F. Liquid Structure at Metal Oxide-Water Interface: Accuracy of a ThreeDimensional RISM Methodology. Chem. Phys. Lett. 2000, 320, 186− 193. (140) Stoyanov, S. R.; Gusarov, S.; Kovalenko, A. Multiscale Modeling of the Adsorption Interaction Between Bitumen Model Compounds and Zeolite Nanoparticles in Gas and Liquid Phase. In Industrial Applications of Molecular Simulations; Meunier, M., Ed.; CRC Press: Boca Raton, FL, 2011; Chapter 14, pp 203−230. (141) Nikolić, D.; Blinov, N.; Wishart, D.; Kovalenko, A. 3D-RISMDock: A New Fragment-Based Drug Design Protocol. J. Chem. Theory Comput. 2012, 8, 3356−3372.

P


Dissipative Particle Dynamics with an Effective Pair Potential from

Recommend Documents