Speeding-up particle simulations of multicomponent polymer systems by coupling to continuum descriptions Marcus Müller1 Institute for Theoretical Physics, Georg-August University, 37077 Göttingen, Germany E-mail: mmueller@theorie.physik.uni-goettingen.de The simulation of structure formation by particle-based simulations poses a computational challenge because of (i) the wide spread of time scales or (ii) the presence of free-energy barriers along the transformation path. A prototypical example of the former difficulty of multiple disparate time scales is the simultaneous presence of stiff bonded interactions, defining the molecular architecture of polymer systems and the weak non-bonded interactions, giving rise to macrophase separation or self-assembly in dense multicomponent systems. A characteristic illustration of the latter problem are nucleation barriers or metastable intermediate states in the course of morphology transformation. Continuum models, in turn, describe the system by a collective order-parameter field, e.g., the composition, rather than particle coordinates, and often do not suffer from these limitations because (i) the stiff molecular degrees of freedom have been integrated out and (ii) advanced numerical techniques, like the string method, exist that identify free-energy barriers and most probable transition paths. Using field-theoretic umbrella sampling, we determine an approximation of the continuum free-energy functional for a specific particle-based model. We illustrate how (i) the on-the-fly string method can identify the minimal free-energy path for the formation of an hourglass-shaped passage (stalk) between two apposing bilayer membranes and (ii) the continuum free-energy functional can be used in conjunction with a heterogeneous multiscale method (HMM) to speed-up the simulation of Lifshitz-Slyozov coarsening in a binary polymer blend by two orders of magnitude. 1 1.1 Soft, coarse-grained particle-based models Length, time, and energy scales in multicomponent polymer melts Soft matter and in particular multicomponent polymer systems are characterized by (i) widely disparate time, length and energy scales, (ii) responsiveness to small driving forces, (iii) a multitude of metastable states, and (iv) structural and chemical complexity of the materials. These challenges require a multiscale approach that often relies on the development and validation of coarse-grained models and the development of new computational strategies. The length, time, and energy scales on the atomic scale, e.g. associated with a covalent bond along the backbone of a polymer, are on the order of Angstrom (bond length), sub-picoseconds (molecular vibrations), and electron Volts (bond energy). The scales associated with a polymer molecule are its root mean-squared end-to-end distance, Re that is on the order of tens of nanometers, the time scale to diffuse its own molecular extension, τ ∼ seconds, and the repulsive interaction free energy between different polymers in a blend, χN kB T ∼ kB T , where kB T denotes the thermal energy scale, χ the Flory-Huggins parameter, and N the number of effective coarse-grained interactions centers along the molecular contour. Length and time scales associated with the collective dynamics of structure formation, i.e. phase separation in a binary homopolymer blend 1 or self-assembly in block copolymer materials, exceed micrometers and hours, respectively. It is quite obvious that no single computational approach can simultaneously address all these different scales and it remains a daunting challenge for a systematic coarsegraining procedure to start out with bond energies of eV and devise a coarse-grained model where free energy differences between effective coarse-grained segments on the order of χkB T ∼ 10−2 kB T ∼ 10−4 eV dictate the qualitative behavior. Additionally, these effective interactions between the coarse-grained segments are free energies, and therefore there is only a limited transferability of the coarse-grained model from one thermodynamic state to another.1 The appropriate choice of the coarse-grained computational models reflects the physical phenomena that one intends to study, i.e. the crystallization of polymers, the glass transition in polymer materials, or phase separation and self-assembly require the coarsegrained description to capture different relevant characteristics of a dense polymer melt. In the following, we will restrict ourselves to structure formation in dense, binary AB polymer materials. These systems are characterized by minute forces that drive structure formation and that cannot yet be adequately predicted by ab initio quantum theory. Therefore the parameters of such models must often be determined directly from experiment. These coarse-grained models describe collective phenomena that can be quantitatively compared to experiments in order to validate the coarse-grained model and, additionally, they provide molecular insights into the structure and dynamics that are often not available experimentally. The wide spread of length, time, and energy scales between the atomistic structure and the morphology imparts a large degree of universality onto the structure formation in multicomponent polymer melts, i.e. systems with different atomistic architectures and interactions exhibit similar behavior on the mesoscopic scale. The appropriate level of description for the study of the mesoscale structure of multicomponent polymer system is the level of an entire macromolecule. On this level of coarse-graining, there are three relevant interactions: (i) bonded interactions, which define the macromolecular architecture, (ii) excluded volume interactions of segments that impart near-incompressibility onto the dense polymer melt, and (iii) repulsion between unlike segment species, which drive the structure formation (i.e. phase separation or self-assembly). These three interactions can be parameterized by three, experimentally accessible, coarse-grained parameters. The length scale of a linear flexible macromolecule, which adopts a Gaussian random-walk configuration, is set is set by Re . The high free-energy costs associated with fluctuations of the total density are set by the inverse isothermal compressibility, κ. Note that in a coarse-grained model it is not necessary to enforce incompressibility down to the scale of a chemical repeat unit or atom but is suffices to limit density fluctuations on the relevant length scale, which is a small fraction of Re . Within √ mean-field theory, the correlation length of density fluctuations is given by ξ = Re / 12κN . The low free-energy scale of interactions between unlike polymer molecules (in a blend) or distinct block in copolymer materials is set by χN in units kB T . This repulsion gives rise to domain formation, and √ the width of the interfaces between domains is given by w0 = Re / 6χN in the limit of large incompatibility. The three coarse-grained parameters, Re , κN , and χN , describe the strengths of the relevant interactions, and they are invariant under changing the number, N , of effective interaction centers that are used to describe the molecular contour. There is one additional, fourth relevant quantity that dictates the behavior of dense 2 
2 polymer melts – the invariant degree of polymerization, N̄ = ρRe 2 /N , where ρ is √ the segment number density. Since in a dense melt, Re = b N , the invariant degree of polymerization is proportional to the number of segments along the chain contour, N̄ = (ρb3 )2 N . The physical meaning of N̄ consists in quantifying the number of neighboring molecules a reference chain interacts with. Since the Gaussian chain conformations are fractal, a Gaussian polymer in three spatial dimensions does not fill space but there are √ on the order N̄ other molecules pervading the volume of the reference chain. This large number of neighbors is one of the important characteristics of dense polymer melts that sets them apart from mixtures of small molecules. In the limit N̄ → ∞, a molecule interacts with infinitely many neighbors and fluctuations of the collective density (or interactions with all the surrounding molecules) are strongly suppressed such that the mean-field theory for polymers – denoted self-consistent field theory – becomes accurate. One important role of computer simulations is to assess the corrections to the mean-field approximation. Likewise, the depth of the correlation hole in the intermolecular pair correlation function, which is important for relating molecular interactions to the Flory-Huggins parameter, or corrections to the Gaussian chain conformations in a dense melt decrease in the limit of large N̄ . Therefore it is important for a coarse-grained model to be able to describe systems with experimentally relevant degree of polymerization. It is important to realize that on this level of coarse-graining one segment corresponds to many chemical repeat units of a chemically realistic representation. While atoms cannot overlap, the centers of mass of a collection of atoms may sit on top of each other. In fact, systematic coarse-graining procedures aiming at reproducing the liquid-like correlations between the coarse-grained segments reveal that the interactions between the coarsegrained segments become the softer the more chemical repeat units a coarse-grained segment represents. As discussed above, the repulsive segmental interactions in the coarsegrained model needs not to be so strong as to enforce incompressibility on the length scale of an atom but we can weaken the repulsive segmental interactions to a level that they are sufficient to suppress density fluctuations on the relevant length scale of a small fraction of Re . This softening of the excluded volume interactions allows for a larger time step in molecular-dynamics simulations or facilitates the use of non-local Monte-Carlo moves (e.g. molecular insertions/deletions via the configuration-bias algorithm). The softness of the interaction is also crucial for representing an experimentally large invariant degree of polymerization, N̄ = 104 . Modeling large values of N̄ = (ρb3 )2 N with particle-based models that include harsh excluded volume interactions between the coarse-grained segments (e.g. lattice models2, 3, 4, 5 or Lennard-Jones potential,6, 7 ) one faces a formidable challenge. The size of a segment, σ, as defined by the range of the harsh √ repulsive interactions, and the statistical segment length of a flexible chain, b ≡ Re / N , are comparable, σ ≈ b. The segment density of a polymer fluid cannot be increased significantly beyond ρσ 3 ≈ 1, because the liquid of segments either crystallizes into a solid or it vitrifies into a glass. Thus, a value of N̄ = 104 requires a large number of segments per 3 2 chain, N = N̄ /(ρb3 )2 ≈ N̄ /(ρσ√ ) ∼ 104 . A small system of linear dimension L = Re 3 is comprised of n = ρL = N N̄ (L/Re )3 ≈ N̄ 3/2 = 106 effective segments. In a dense melt, these long entangled chains will reptate,8, 9 and the time to diffuse a distance Re scales like τ = τ0 N 3 ∼ N̄ 3 where τ0 is a N -independent microscopic time scale. To follow the system over one characteristic time one needs about N̄ 9/2 = 1018 elementary 3 segment moves. Contrary, if the harsh excluded volume interaction is replaced by a soft repulsion, one will eliminate the constraint ρo b3 . 1, because solidification or vitrification √ can be avoided. In this case, one can choose a much larger segment density, ρb3 ∼ N̄ . For instance, choosing ρb3 = 18, we can model a value of N̄ = 104 by using N = 31 segments along the molecular contour. This discretization of the molecular architecture is still sufficient to capture the characteristics of the random-walk-like conformations on the scale Re . Within the soft coarse-grained model, a system of size L = Re contains only 3 200 segments. Moreover, these non-entangled polymers obey Rouse dynamics with a re√ laxation time τ = τo N 2 . Thus the simulations require only N 3 N̄ ≈ 3 · 108 elementary moves, which is 11 orders of magnitude less than in coarse-grained models, where excluded volume is enforced on the scale of a segment. For this reason, soft coarse-grained models are very efficient in describing polymer systems with a realistically large value of N̄ and allow us to study collective phenomena on the length scale of Re and beyond. This ability can be traced back to the rather coarse representation of the molecular contour and the concomitant large number of monomeric repeat units that are lumped into an effective coarse-grained segment. In order to identify the length and time scales of the soft coarse-grained model let us consider a melt of polystyrene with a molecular weight of Mw = 100 000 or 962 chemical repeat units C8 H8 . The statistical segment length of a chemical repeat unit is about 0.7 nm yielding Re = 21.7 nm. A mass density of 1.06 g/cm3 translates into an invariant degree of polymerization of N̄ ≈ 4200. Using a typical self-diffusion coefficient of D = 10 nm2 /s,10 we obtain a characteristic time scale of τ = Re 2 /D = 47 s. In computer simulations of soft, coarse-grained models one can study systems with a few million coarse-grained segments. Assuming a chain discretization of N = 32, i.e. one coarse-grained segment correspond to 30 chemical repeat units, a typical system is comprised of some 30 000 molecules corresponding to a linear extension L ∼ 8Re ∼ 0.17 µm of a cubic system. A typical simulation with 106 elementary steps per segment corresponds to 100 τ or 1.5 hours. Thus soft, coarse-grained models are able to reach the experimentally relevant length and time scales of phase separation and self-assembly in polymer blends and block copolymer materials. 1.2 Soft, coarse-grained particle-based models for multicomponent polymer melts We use a minimal, soft, coarse-grained model that captures the three relevant interactions. In the following, we distinguish between bonded interactions, which define the molecular shape and its fluctuations, and non-bonded interactions, that impart near-incompressibility onto the dense melt and drive the structure formation. Since a coarse-grained segments is comprised of many chemical repeat units, the distance distribution between neighboring coarse-grained segments along the macromolecule is Gaussian due to the central limit theorem and orientational correlations along the backbone of the chemical structure have decayed on the length scale of a coarse-grained segment. Therefore, we model the universal aspects of the molecular shape by the discretized Edwards Hamiltonian. N −1 X 3(N − 1) Hb [ri (s)] 2 = [ri (s + 1) − ri (s)] , 2 kB T 2R e s=1 4 (1) where we consider n polymers with N segments in a volume V . {ri,s } with i = 1, · · · , n and s = 1, · · · , N denotes the set of segment coordinates that completely specifies the configuration of our system. The density of the melt is ρ = nN/V and Re denotes the root mean-squared end-to-end distance of an ideal chain, i.e. in the absence of non-bonded interactions. The soft, pairwise interactions can be re-written in the form of a free-energy functional11 Hnb ({r}) = Fnb [φ̂A (r|{r}), φ̂B (r|{r})] (2) where the local microscopic densities, φ̂A and φ̂B , depend on the particle coordinates, {r}. 1 X δ(r − riA ) (3) φ̂A (r|{r}) = ρ0 i A The sum runs over all A segments irrespectively to which molecule they belong. A typical local free-energy functional for non-bonded interactions in an AB binary melt can be written as   Z dr χN κN Fnb [φ̂A , φ̂B ] p 2 2 = N̄ ( φ̂ + φ̂ − 1) − ( φ̂ − φ̂ ) (4) A B A B kB T 2 4 Re 3 where χ is the bare Flory-Huggins parameter, and κ is the bare, dimensionless, inverse isothermal compressibility. Like the end-to-end distance, the actual energy of mixing or compressibility slightly deviates from the parameters of the Hamiltonian due to fluctuation/correlation effects. The advantage of this formulation is that it offers a general strategy to systematically incorporate thermodynamic information into the soft, coarse-grained model. Eqs. (3) and (4) are not suitable for computer simulation; the δ-function needs to be regularized. Either one computes the local densities by smearing the δ-function out over a volume ∆L3 or one employs a collocation lattice of√grid spacing ∆L. Typically, ∆L is comparable to the statistical segment length, b = Re / N of the coarse-grained model and smaller than the width of the AB interfaces, w0 . In the first method, one represents the δ-function in Eq. (3) as a limit of a weighting function, ω, and defines a weighted density12 Z 3 0 d r 1 X φ̂A,ω (r|{r}) = ω(|r − r0 |)φ̂A (r0 |{r}) = ω(|r − riA |) (5) 3 ∆L ρ∆L3 i A R with normalization d3 r ω(|r|) = ∆L3 . In the simplest case, ω is proportional to the characteristic function of a sphere. A quadratic term in the excess free energy yields a density-dependent pairwise potential.13, 14 p Z d3 r 1 X N̄ v(riA , rjB ) (6) 3 φ̂A (r|{r})φ̂B (r|{r}) = N 2 Re i ,j A B R d3 r Re 3 which is translationally invariant and isotropic, i.e., v(|r0 − r00 |) = √1 ∆L 3 ∆L3 ω(|r − N̄ R d3 r 0 00 r |)ω(|r − r |). The N̄ -dependence of the integrated strength, Re 3 v(|r|) = √1 , guarN̄ antees that the limit of high density remains well defined. 5 In the grid-based scheme, one discretizes space in cubic cells of linear dimension, ∆L. Each cell is identified by its index, c. We define the local microscopic densities on the grid by assigning particle positions to the grid cells according to15, 13 1 X d3 r Π(c, r) φ̂ (r) = Π(c, ri,s ) (7) A ∆L3 ρ∆L3 i A R P The assignment function fulfills c Π(c, r) = 1 ∀r and d3 r Π(c, r) = ∆L3 ∀c i.e. the contribution of a particle to all cells adds up to unity irrespectively of its position, and the volume assigned to each grid cell is ∆L3 . In the simplest case, Π(c, r) is the characteristic function of a grid cell. The grid-based method also yields pairwise interactions according to Eq. (6) but, since they make reference to the underlying lattice, they are no longer translationally and rotationally invariant, v(r0 , r00 ) = 3 P 0 00 √1 Re 3 c Π(c, r )Π(c, r ). Therefore, one needs to resort to special simulation techN̄ ∆L niques for computing the pressure and care hat to be exerted to control the effect of selfinteractions.16 However, in the grid-based approach, the energy of a segment with its surrounding can be very efficiently computed from the knowledge of the density on the collocation lattice. In the former weighting-function method, in turn, the energy involves the explicit computation of the pairwise interactions between a segment and its neighbors. This calculation is performed via a cell list, where the cell’s linear dimension is the range of the pairwise interaction, O(∆L). All interactions in the 27 cells around the one that contains the segment have to be considered. For a typical choice of parameters, N = 32, N̄ = 104 , ∆L/Re = 1/6 this amounts to O(102 ) interaction pairs. Thus the grid-based technique offers a significant computational advantage for dense polymer systems. Z φ̂A (c|{r}) = 1.3 Strong bonded and weak non-bonded forces and SCMF Simulations Due to the computational speed-up we use the grid-based approach in the following. Since the pairwise interactions are not translationally invariant, we explore the configuration space of the soft, coarse-gained model by Monte-Carlo simulations. It is worth noting that for fine discretization of the molecular contour, N  1, there is a pronounced difference between the strong bonded forces, fb , that define the molecular architecture and the weak non-bonded forces, fnb , that drive structure formation. kB T kB T √ χkB T 1 kB T p fb ∼ ∼ and fnb ∼ (8) · N ∼ · 6(χN )3 · b Re w0 Re N i.e. fb /fnb ∼ N 3/2 . In molecular dynamics simulations, one would use multiple time-step integrators (rRESPA)17 to cope with this disparity of forces. In Monte-Carlo simulations, one can use the Single-Chain-in-Mean-Field (SCMF) algorithm15, 18 to exploit the separation between the strong, rapidly fluctuating, bonded interactions, which dictate the size of a segmental movement in one Monte Carlo step, and the weak, non-bonded interactions, which only very slowly evolve in time. In SCMF simulations, we temporarily replace the pairwise interactions, Eq. (2), of a segment with its HSCMF nb surroundings by the interaction of a segment with an external field, i.e. = kB T h i 3 P ρ∆L c wA (c)φ̂A (c|{r}) + wB (c)φ̂B (c|{r}) , where the external field, wA /N that N 6 acts on A segments is frequently calculated from the local fluctuating densities according ∂Fnb N to wA (c) = ρ∆L 3 ∂φ (c) . A SCMF simulation cycle is comprised of two parts: 1) evolve A the polymer conformations in the external fields, wA and wB , for a small, fixed amount of Monte-Carlo steps. We employ Smart-Monte-Carlo moves, using the strong bonded forces to bias the proposal of a trail displacement.19, 19 During this Monte-Carlo simulations the molecules do not interact with each other and the simulation of independent chain molecules can be straightforwardly implemented on parallel computers. 2) recalculate the external fields from the instantaneous densities. In this second step, fluctuations and correlations are partially restored. Then the simulation cycle commences again. The quasiinstantaneous field approximation that consists in replacing the interactions via frequently updated, fluctuating, external fields will be accurate, if the change of the local composition between successive updates of the external fields is small. This property is controlled 1 by the parameter, ε = N ρ∆L 3 , which plays a similar role as the Ginzburg parameter in a mean-field calculation. In contrast to the Ginzburg parameter, however, ε depends on the discretization of space, ∆L, and molecular contour, N , and these parameters are chosen such that the quasi-instantaneous field approximation is accurate.13 1.4 Barrier and time-scale problem of particle-based models In spite of the benefits of soft, coarse-grained models, the investigation of the kinetics of phase separation or self-assembly in computer simulations of particle-based models is computationally demanding. Two effects contribute to this difficulty: (i) barrier problem – In the course of structure formation, multiple free-energy barriers must be overcome. Since collective structure formation involves many molecules, free-energy barriers typically exceed kB T , and rare thermal fluctuations are required to overcome them. For the favorable case in which it is possible to identify a suitable and simple reaction coordinate, or when one can identify a low-dimensional subspace that characterizes the barriers, a variety of computational techniques have been devised to “flatten” the free-energy landscape and facilitate the exploration of phase space or to compute the saddle-points of the free-energy landscape that dictate the kinetics of structure formation.20, 21 (ii) time-scale problem – Even if the time evolution is completely down hill in free energy, the kinetics of the order parameter can be intrinsically slow because the thermodynamic driving force does not efficiently generate a concomitant current. A prototypical example is the diffusion of one molecule from one domain to another, as it occurs in Lifshitz-Slyozov coarsening in binary blends,22 the diffusion across lamellae in symmetric block copolymers, or the exchange of amphiphiles between micellar aggregates. In these cases, molecules have to “tunnel” through an unfavorable domain, and this thermally activated process dramatically slows down the current. These two types of problems can be addressed by concurrent coupling of the particlebased model to a continuum description. 7 2 2.1 Continuum models Order parameter and free-energy landscape In a continuum approach, the system configuration is entirely described through a collective order-parameter, i.e. a continuum field that does not make references to the properties of individual molecules. The choice of the order parameter is critical and crucially depends on the problem at hand. In the following we consider examples where the difference between the local A and B densities provides an appropriate order parameter, m(r) = φA (r) − φB (r). Then, one can define the free energy as a functional of the order parameter m(r) via the trace over all particle conformations compatible with m(r) Z QnAB QN h i H({r}) Y F [m] −k T i=1 t=1 dri,t − kB T B ≡ e δ m(r) − φ̂ (r|{r}) + φ̂ (r|{r}) (9) e A B AB N nAB !λ3n T r where we considered nAB molecules consisting of N segments. λT is the thermal deBroglie wavelength, and H({r}) denotes the interactions of the underlying particle-based model. Eq. (9) guarantees that the partition functions of the particle-based modelh and of ithe QN H({r}) R R Qni=1 AB − t=1 dri,t continuum description are identical e kB T = Dm exp − FkB[m] 3nAB N T . nAB !λT Knowledge of the free-energy functional (or landscape) allows one to draw important conclusions: (i) Within the mean-field approximation, minima of F [m] correspond to (meta)stable states. (ii) If there is a clear separation of time scales between the fast singlechain dynamics and the slow kinetics of the order-parameter, the molecular conformations will be in equilibrium with the instantaneous order-parameter, i.e., they sample the equilibrium distribution that is compatible with the order-parameter field, Eq, (9). In this limit, the qualitative features of the order-parameter dynamics can be inferred from the freeenergy landscape. Most importantly neglecting thermal fluctuations, one can distinguish two types of structure formation kinetics – spinodal self-assembly or phase separation and nucleation. Having identified the free-energy landscape as a function(al), F [m], of a slowly evolvδF ing order parameter, one can compute the thermodynamic force, ∇ δm(r) , that drives structure formation. The Onsager coefficient, Λ, connects this thermodynamic force to the current of the order parameter: Z δF (10) j(r) = − dr0 Λ(r, r0 )∇0 δm(r0 ) Since the order parameter is related to the densities of the different species, it is conserved and obeys the continuity equation23, 24, 25 ∂m = −∇j (11) ∂t In an incompressible systems, the currents of A and B segments cancel and Λ ∼ φA φB = (1 − m2 )/4. It is this factor that gives rise to intrinsically slow dynamics, cf. Sec. 1.4, i.e. in strongly segregated systems the kinetics can be protracted even if there is a strong thermodynamic driving force. a a General expressions for relating the Onsager coefficient to the dynamics of the underlying macromolecules 8 Eqs. (10) and (11) can be augmented by random noise terms such that the dynamics is able to overcome barriers in the free-energy landscape.27, 28 Thermal fluctuations are often neglected and the mean-field approximation is invoked; in order to address fluctuation effects one has to cope with short-length scale fluctuations, which lead to UV-divergencies. Different simple forms of the free-energy functional, F [m], have been proposed on the basis of general symmetry principles. A common description of binary blends is provided by the Ginzburg-Landau square-gradient functional.29 Microphase separation of block copolymer materials can be described by the Otha-Kawasaki functional30 or the SwiftHohenberg approach.31 The small number of parameters that enter such a continuum description can be qualitatively related to physically accessible quantities like the segregation inside the domains or the intrinsic widths of interfaces. Because they ignore all molecular degrees of freedom, these continuum models are computationally efficient. Additionally, sophisticated methods have been devised to identify barriers and minimal free-energy paths32 , and the effects of small Onsager coefficients can be mitigated by using a large time step for integrating Eq. (11). 3 Systematic parameterization of a continuum model: field-theoretic umbrella sampling and force matching The barrier and time-scale problem in particle-based models can be addressed by coupling them to a continuum model in the framework of the heterogeneous multiscale method (HMM).33, 34 To this end, one has to estimate the free-energy functional, F [m], of the particle-based model. Two computational strategies have been devised to this end: fieldtheoretic umbrella sampling35 and field-theoretic force matching.16 In both cases, one does not directly obtain the free-energy functional but rather the chemical potential, µ(r|m) ≡ δF δm(r) , for a specific configuration of the order parameter. In field-theoretic umbrella sampling,35 one adds to the interactions of the particle-based model an umbrella potential that restrains the local microscopic densities, φ̂A (r|{r}) − φ̂B (r|{r}), to the local value of the order-parameter, m(r), at each point in space Z  i2 λh Hfup ({r}) = dr m(r) − φ̂A (r|{r}) − φ̂B (r|{r}) (12) 2 The integral in Eq. (12) is evaluated using a collocation lattice (see. Sec. 1.2). In the limit, λ → ∞, the Boltzmann factor of this field-theoretic umbrella potential converges to the δ-function constraint in Eq. (9) that projects out the microscopic particle configurations compatible with the order parameter,35, 36 and free energy of the restrained system with the field-theoretic umbrella potential, Fλ [m], converges towards the constraint free-energy. The chemical potential can be calculated according to have been devised.26 For the Rouse model with inverse friction * Λ(r, r0 ) ≈ 1 ζ = ND kB T one obtains + X ∂ 1 ∂ ND g(r, r0 ) [φ̂A (r|{r}) − φ̂B (r|{r})] · · [φ̂A (r0 |{r}) − φ̂B (r0 |{r})] ≈ (1−m2 ) ∂r ζ ∂r ρk T V i i B i where the last expression refers to a symmetric homopolymer blend in the disordered state and the non-locality is characterized by the single-chain correlation function, g(r). 9 δFλ = δm(r)  δHfup δm(r)  δF = µ(r|m) δm(r) λ (13) where the average h· · · iλ is performed in the restrained system. Independent from λ, this average has to be sampled for about one molecular relaxation time, τ to accurately calculate the local chemical potential.35 In field-theoretic force matching, one alternatively can use the thermodynamic relation between the force, KA (r), acting on A-segments in a volume element around position, r, and the gradient of the excess chemical potential:16 µλ (r|m) ≡ D  E = λ m − φ̂A − φ̂B λ→∞ → λ ρ hKA (r)i|m(r) = −∇ [µA (r|m) − ρkB T ln φA (r)] (14) and µ = µA − µB . The force is determined in configurations that are characterized by the order parameter, m(r). The advantage of this technique is that it does not rely on the limit λ → ∞ or the use of a collocation grid. For polymers, however, there are large cancellation effects of forces similar to the atomistic expressions for the virial pressure. 4 4.1 Applications Barrier problem: Minimum free-energy path (MFEP) of stalk formation In order to overcome the barrier problem and find a suitable path along which structure formation proceeds, one can adopt an equation-free approach, where no Ansatz for the explicit form of the free-energy functional is required.37 Knowing the chemical potential µ(r|m), we use the string method32 to find the minimal free energy path (MFEP) that connects the starting and ending order-parameter configurations.38, 39, 40 The MFEP is a string of morphologies, ms (r), where s denotes the contour parameter along the string and the squared distance between R two neighboring morphologies, ms (r) and ms0 (r), along the string is given by ∆2s,s0 ∝ dr [ms (r) − ms0 (r)]2 . The MFEP is defined by the condition that the thermodynamic force in the direction perpendicular to the path vanishes R dms (r0 ) dms (r) dr0 µ(r0 |ms ) ds (15) ∇⊥ F[ms ] = µ(r|ms ) −  2 R 0 ds s (r ) dr0 dmds Thus the defining condition for the MFEP can be solely expressed by the chemical potential that we obtain in the particle-based model via field-theoretic umbrella sampling, Eq. (13). The MFEP is efficiently determined numerically by the improved string method,32, 38 which consists of a two-step cycle: (i) F is minimized by evolving the morphologies according [m] to ∆ms (r) = −µ(r|ms )∆ε with µ(r|m) = δF δm(r) ; and (ii) ms (r) is re-parameterized via a third-order spline at each point, r, to restore uniform distance of the morphologies along the string. One application is illustrated in Fig. 1, where this techniques has been employed to study the formation of an hour-glass shaped, hydrophobic passage (stalk) between two apposing lamellar sheets in a copolymer-homopolymer mixture.40 By virtue of the universality of the structure of amphiphilic systems,41 this model can be conceived as a representation of lipid membranes – the A and B blocks corresponding to hydrophobic tails and hydrophilic heads of lipid molecules and the B-homopolymers representing the solvent. 10 Figure 1. Minimum free-energy path (MFEP) obtained by the on-the-fly string method of stalk formation. Adapted from Ref.40 The MFEP, msi (r), is discretized into 24 particle-based systems and intermediate values of s are obtained by point-wise spline interpolation. The free energy along the MFEP R [ms ] s (r) δF [ms ] ∗ is obtained by dFds = dr ∂m∂s δm(r) . and the transition state, m , is identified as [ms ] the maximum on the MFEP, dFds = 0. The left axis of Fig. 1 presents the free energy, F [ms ], along the MFEP in units of γd2 , √ p 2 where γRe / N̄ kB T ≈ χ0 N/6 and d ≈ 1.82Re denote the AB (oil-water) interface tension and the lamella (bilayer) thickness, respectively. Typical experimental values of lipid membranes are d ≈ 3.6nm and γd2 = 155kB T . The contour plots depict cross sections of the order parameter, ms (r), for the stable, apposing-bilayer morphology and the metastable stalk morphology. The snapshot depicts a particle configuration restrained by the field-theoretic umbrella potential, Eq. (12), using the order parameter, ms∗ (r), at the saddle point, s∗ = 0.532, of the MFEP. Hydrophilic beads are colored yellow, hydrophobic beads are shown in red, solvent (homopolymer) particles are not shown. Only every 10th copolymer is depicted corresponding to a typical density in a lipid system. This static information is complemented by the probability that configurations along the MFEP have transformed in the course of simulations into two apposed bilayers at a specified time after the restraining field-theoretic umbrella potential has been removed (blue, right axis). Results have been obtained for 256 independent configurations at each 11 Figure 2. Illustration of one cycle of HMM for the initial stages (spinodal decomposition) of phase separation in a homopolymer blend. Adapted from Ref.43 value of s. The probability is a sharply varying function along the MFEP and the position, s, at which there is a 50 − 50 chance of reaching either the morphology of two apposing bilayers or the stalk, agrees with the saddle-point of the MFEP. The hydrophobic bridge that connects the two lamellae – denoted by stalk – has attracted much interest in the context of bilayer membrane fusion.42 Our particle-based simulations provide direct microscopic insights into the transition state that consists of only few hydrophobic tails that bridge between the bilayers and that constitutes a free-energy barrier of 16kB T in a lipid system. 4.2 Time-scale problem: Heterogeneous multiscale method applied to Lifshitz-Slyozov coarsening in a binary polymer blend By knowing the free-energy functional, F [m], one can also mitigate the time-scale problem by concurrently coupling the particle model to the corresponding continuum model. This HMM33, 34, 35 comprises three steps, which are illustrated in Fig. 2: 1) estimate the parameters of the continuum description, 2) propagate the continuum model for a large time step, ∆t, and 3) seamlessly generate a new particle configuration compatible with the new order-parameter field. Then the cycle commences again. In step 1), one needs to compute the chemical potential, µ, and the Onsager coeffi- 12 cient, Λ, for a specific configuration of the particle model. The former can be obtained by field-theoretic umbrella sampling or field-theoretic force matching. This equation-free strategy, however, would require the chemical potential be frequently computed because the local chemical potential significantly changes on the time scale where an AB interface has moved a distance comparable to its intrinsic width. Therefore, rather than using the chemical potential directly, it is useful to exploit the knowledge of µ(r|m) to parameterize an explicit Ansatz, Ftrial [m], for the free-energy landscape that contains a small number of variational parameters, {α}, like e.g. the Flory-Huggins parameter and the coefficient of the square-gradient term. Using the measured µ, one can adjust these parameters to minimize the deviation between F and Ftrial , i.e. we choose {α} such that 2  R trial [m] → min. dr µ(r|m) − δFδm(r) This Ansatz tacitly assumes that Ftrial [m] with the same set of {α} is able to describe the entire system, e.g. it can simultaneously describe the composition inside large domains and the profiles across AB interfaces. If the Ansatz were perfect, the parameters, {α}, would not depend on the order-parameter configuration, and one could employ the onceparameterized free-energy functional to predict the entire kinetics of structure formation of the particle-based model. In practice, the optimal parameters will slightly depend on the specific m(r). For instance, we anticipate changes of Ftrial when the segregation of the domains changes or the structure of the AB interfaces is altered. The residual minimum indicates the quality of the Ansatz, Ftrial , signals the need for re-parameterization, and allows for a systematic improvement of Ftrial by including additional terms. Moreover, the computational time required for computing the small number of parameters, {α}, is significantly smaller than accurately computing the chemical potential at each point in space because one can substitute the time average of a local quantity by a spatial average over the entire system. It is important to realize that changes of the thermodynamic state that require reparameterization occur on a time scale that is much longer than the motion of interfaces. Hence, Ftrial [m] can predict the structure formation for a much larger time interval than µ(r|m), and the time step, ∆t, of a single cycle can be significantly larger in HMM than in an equation-free scheme that directly uses the local chemical potential.44 Since the continuum model is not explicitly concerned with the stiff bonded degrees of freedom the time scale can be adjusted to the intrinsically slow process and step 2) of the HMM scheme takes a vanishingly small computation time compared to the propagation of the particle-based model. To seamlessly generate a new particle configuration in step 3), which corresponds to the new order-parameter field, m(r, t+∆t), we use the same field-theoretic umbrella potential that has been employed to compute the chemical potential. Using the new order-parameter at time t + ∆t in the field-theoretic umbrella potential, Eq. (12), one creates a large thermodynamic force, −λkB T ∇[m(r, t + ∆t) − (φ̂A − φ̂B )] towards the new order-parameter configuration. This strong force amplifies the weak thermodynamic driving forces of the non-bonded interactions in the original model by a factor that is proportional to the strength λ, speeding-up the relaxation towards m(r, t + ∆t) compared to the original dynamics of structure formation. This rational suggests that λ be chosen as large as possible in order to achieve the maximal speed-up. There are, however, two limitations: (i) The thermodynamic force of the field-theoretic umbrella potential should amplify the weak thermody- 13 namic driving force of the original model but they must not exceed the strong non-bonded forces that dictate the single-molecule dynamics. Otherwise, the underlying particle dynamics will be altered and the time step used to evolve the particle-based model has to be reduced. (ii) The estimate of the speed-up relies on linear response theory that fails already at moderately large values of λ. In this case, one might need an additional time τ to relax the molecular conformations to the equilibrium statics within the field-theoretic umbrella potential. This step 3) also provides a strategy for computing the Onsager coefficient from the last stage of relaxation towards the new equilibrium in the restrained system. In the case of large λ the field-theoretic umbrella potential dominates the thermodynamic force and, since the difference m − (φ̂A − φ̂B ) is small, linear response theory is appropriate and predicts an exponential relaxation towards the restrained equilibrium. Alternatively, one can estimate Λ by comparing the kinetics of structure formation of the original particlebased model with the prediction of continuum approach. One can additionally speed-up the relaxation towards the new order-parameter field by computing the average current, j̄, during the time interval, ∆t, from the continuum model and estimate the concomitant time-averaged velocity fields, v̄A (r) and v̄B (r). Then, one couples these flow fields to the particle model via an additional drag force, Fi = γ v̄A (ri ) acting on an A particle at position, ri . γ is a friction coefficient. Applying this force at the initial stage of relaxation towards the new order parameter, one accelerates the generation of a new particle configuration. Also in this case, an additional molecular relaxation time without flow is required to bring the molecular conformations into equilibrium with the field-theoretic umbrella potential. Steps 1) and 3) of HMM require a time of the order of the molecular relaxation time, τ . The computational cost of propagating the continuum model is negligible and thus the computational speed-up is of the order ∆t/τ . Using an accurate free-energy functional, Ftrial , that is suitable for describing the slow structure formation over a long time interval, ∆t, without the need for re-parameterization, large speed-ups are feasible. The sogenerated particle configurations can subsequently be used to investigate the single-chain conformations and dynamics, which is not accessible in the continuum model. One application of HMM is illustrated in Fig. 3 where the evaporation of chains from a drop in a binary AB homopolymer blend is investigated. The system of geometry 12Re × 6Re ×6Re is comprised of two A domains – a spherical drop with excess ∆A of A segments and a planar slab-like domain that spans the system via the periodic boundary conditions. The morphologies are illustrated in the inset images. Due to the curvature of the drop’s interface, the chemical potential is inside the drop is higher than in the planar domain and A molecules evaporate from the drop and condense onto the planar domain (LifshitzSlyozov coarsening22 ). This process is protracted because the Onsager coefficient inside the strongly segregated B-rich matrix is very small. Fig. 3 depicts the linear shrinking of the drop’s volume with time. The red solid line corresponds to the simulation of the particle-based model and the black dashed line depicts the prediction of the continuum model. The continuum model is a Ginzburg-Landau square-gradient model where we have adjusted the effective incompatibility and the coefficient of the square-gradient term. The Onsager coefficient was determined by comparing the time evolution of the particle-based model and the continuum model at early times. The so-parameterized continuum model accurately describes the entire drop evaporation. The steep doted lines show the relaxation 14 Figure 3. HMM of macrophase separation in a soft, coarse-grained model of an AB homopolymer blend, χN = 5. The inset presents an enlargement of the main panel. Snapshots illustrate configurations of the particle-based model (right) at t0 + 41τ and t0 + 82τ along the original time evolution and the corresponding configurations after the relaxation step 3) of one HMM-cycle. (left). Adapted from Ref.44 of the particle-based model towards the new order parameter at a later time ∆t using the field-theoretic umbrella potential and an initial coupling to j. The green and black lines with symbols present the free time evolution of the particle-based model restarted after one HMM step. The unrestrained particle simulation restarted with the new configuration captures the behavior of the original trajectory after time ∆t indicating that the HMM scheme also captured the decay of the composition across the B-matrix, which dictates the evaporation rate. In this example speed-ups of ∆t/τ = 41 and 82 are achieved with respect to SCMF simulations of the particle-based model.44 5 Concluding Remarks Soft, coarse-grained models are well suited to efficiently investigate the universal equilibrium behavior of multicomponent polymer blends and copolymer materials in the liquid state. They can successfully address the relevant time and length scales of structure formation and allows us to systematically explore the structural and chemical diversity of multicomponent materials and provide structural and dynamic insights on the molecular level that are often not readily available in experiments. These models are a good start- 15 ing point for investigating the collective dynamics of phase separation and self-assembly in nanostructured materials. Given the multitude of metastable states, there is a great potential in controlling and directing the dynamics of structure formation and identifying mechanisms of collective structure transformations. Due to the widely spread time and length scales, understanding and reliably predicting the practically important relation between the single-molecule dynamics and the kinetics of morphological changes remains a formidable challenge and computational techniques that seamlessly couple different levels of description will be instrumental in exploring how the collective dynamics can be tailored by the underlying motion of the molecules or application of external fields. Acknowledgments I am indebted to K.Ch. Daoulas, A.-C. Shi, J.J. de Pablo, Y. Smirnova, S. Glatz, G. Marelli, and M. Fuhrmans for fruitful and enjoyable collaborations. Financial support by the Volkswagen foundation and the SFB 803-B3 are gratefully acknowledged. I thank the Jülich Supercomputer Center and the HLRN Hannover/Berlin for generous allocation of computational resources. References 1. A. A. Louis, Beware of Density Dependent Pair Potentials, J. Phys.: Condens. Matter, 14, 9187–9206, 2002. 2. P. J. Flory, Thermodynamics of High Polymer Solutions, J. Chem. Phys., 9, 660, 1941. 3. M. L. Huggins, Solutions of Long Chain Compounds, J. Chem. Phys., 9, 440, 1941. 4. I. Carmesin and K. Kremer, The bond fluctuation method: a new effective algorithm for the dynamics of polymers in all spatial dimensions, Macromolecules, 21, 2819– 2823, 1988. 5. M. Müller, Miscibility Behavior and Single Chain Properties in Polymer Blends: a Bond Fluctuation Model Study, Macromol. Theory Simul., 8, 343–374, 1999. 6. G. S. Grest and K. Kremer, Molecular dynamics simulations for polymers in the presence of a heat bath, Phys. Rev. A, 33, 3628, 1986. 7. M. Murat, G. S. Grest, and K. Kremer, Statics and Dynamics of Symmetric Diblock Copolymers: a Molecular Dynamics Study, Macromolecules, 32, 595–609, 1999. 8. P. G. de Gennes, Reptation of a polymer chain in presence of fixed obstacles, J. Chem. Phys., 55, 572, 1971. 9. M. Doi and S.F. Edwards, The Theory of Polymer Dynamics, Oxford University Press, New York, 1994. 10. Markus Antonietti, Jochen Coutandin, and Hans Sillescu, Chain length and temperature dependence of self-diffusion coefficients in polystyrene, Makromol. Chem., Rapid Comm., 5, 525–528, 1984. 11. K. C. Daoulas and M. Müller, Comparison of Simulations of Lipid Membranes with Membranes of Block Copolymers, Adv. Polym. Sci., 224, 197–233, 2010. 12. M. Laradji, H. Guo, and M. J. Zuckermann, Off-lattice Monte Carlo simulation of polymer brushes in good solvent, Phys. Rev. E, 49, 3199, 1994. 16 13. K. Ch. Daoulas and M. Müller, Single Chain in Mean Field simulations: Quasiinstantaneous field approximation and quantitative comparison with Monte Carlo simulations, J. Chem. Phys., 125, 184904, 2006. 14. F. A. Detcheverry, D. Q. Pike, P. F. Nealey, M. Müller, and J. J. de Pablo, Monte Carlo simulation of coarse grain polymeric systems, Phys. Rev. Lett., 102, 197801, 2009. 15. M. Müller and G. D. Smith, Phase separation in binary mixtures containing polymers: a quantitative comparison of Single-Chain-in-Mean-Field simulations and computer simulations of the corresponding multichain systems, J. Polym. Sci. B: Polymer Physics, 43, 934–958, 2005. 16. M. Müller, Studying amphiphilic self-assembly with soft coarse-grained models, J. Stat. Phys., 145, 967, 2011. 17. M. Tuckerman, B.J. Berne, and G.J. Martyna, Reversible multiple time scale molecular dynamics, J. Chem. Phys., 97, 1990, 1992. 18. K. Ch. Daoulas, M. Müller, J. J. de Pablo, P. F. Nealey, and G. D. Smith, Morphology of Multi-Component Polymer Systems: Single-Chain-in-Mean-Field Simulation Studies, Soft Matter, 2, 573–583, 2006. 19. P. J. Rossky, J. D. Doll, and H. L. Friedman, Brownian Dynamics as Smart MonteCarlo Simulation, J. Chem. Phys., 69, 4628–4633, 1978. 20. F. G. Wang and D. P. Landau, Efficient multiple-range random walk algorithm to calculate the density of states, Phys. Rev. Lett., 86, 2050–2053, 2001. 21. Alessandro Laio and Michele Parrinello, Escaping free-energy minima, Proc. Natl. Acad. Sci. USA, 99, 12562–12566, 2002. 22. I. M. Lifshitz and V. V. Slyozov, The kinetics of precipitation from supersaturated solid solutions, J. Phys. Chem. Solids, 19, 35, 1962. 23. J. W. Cahn, Phase separation by spinodal decomposition in isotropic systems, J. Chem. Phys., 42, 93–99, 1965. 24. P. C. Hohenberg and B. I. Halperin, Theory of dynamic critical phenomena, Rev. Mod. Phys., 49, 435–479, 1977. 25. U. M. B. Marconi and P. Tarazona, Dynamic Density Functional Theory of Fluids, J. Chem. Phys., 110, 8032–8044, 1999. 26. Kyozi Kawasaki and Ken Sekimoto, Dynamical theory of polymer melt morphology, Physica A, 143, no. 3, 349 – 413, 1987. 27. H.E Cook, Brownian motion in spinodal decomposition, Acta Metallurgica, 18, 297 – 306, 1970. 28. R. Petschek and H. Metiu, Computer simulation of the time-dependent GinzburgLandau model for spinodal decomposition, J. Chem. Phys., 79, 3443–3456, 1983. 29. J. W. Cahn and J. E. Hilliard, Free energy of a nonuniform system I: Interfacial free energy, J. Chem. Phys., 28, 258–267, 1958. 30. T. Ohta and K. Kawasaki, Equilibrium Morphology of Block Copolymer Melts, Macromolecules, 19, 2621–2632, 1986. 31. J. Swift and P. C. Hohenberg, Hydrodynamic fluctuations at the convective instability, Phys. Rev. A, 15, 319–328, 1977. 32. W. E, W. Ren, and E. Vanden-Eijnden, Simplified and improved string method for computing the minimum energy paths in barrier-crossing events, J. Chem. Phys., 126, 164103, 2007. 33. W. E, B. Engquist, X. T. Li, W. Q. Ren, and E. Vanden-Eijnden, Heterogeneous Mul- 17 tiscale Methods: A Review, Comm. Computational Phys., 2, 367–450, 2007. 34. W. E., W. Q. Ren, and E. Vanden-Eijnden, A general strategy for designing seamless multiscale methods, J. Comp. Phys., 228, 5437–5453, 2009. 35. M. Müller, Concurrent coupling between a particle simulation and a continuum description, Eur. Phys. J. Special Topics, 177, 149, 2009. 36. L. Maragliano and E. Vanden-Eijnden, A temperature accelerated method for sampling free energy and determining reaction pathways in rare events simulations, Chem. Phys. Lett., 426, 168–175, 2006. 37. I. G. Kevrekidis, C. W. Gear, and G. Hummer, Equation-Free: The computer-aided analysis of complex multiscale systems, AICHE J., 50, 1346–1355, 2004. 38. L. Maragliano and E. Vanden-Eijnden, On-the-Fly string method for minimum free energy paths calculation, Chem. Phys. Lett., 446, 182–190, 2007. 39. Magdalena Venturoli, Eric Vanden-Eijnden, and Giovanni Ciccotti, Kinetics of phase transitions in two dimensional Ising models studied with the string method, J. Math. Chem., 45, 188–222, 2009. 40. M. Müller, Y. G. Smirnova, G. Marelli, M. Fuhrmans, and A. C. Shi, Transition path from two apposed membranes to a stalk obtained by a combination of particle simulations and string method, Phys. Rev. Lett., 108, 228103, 2012. 41. M. Müller, K. Katsov, and M. Schick, Coarse-Grained Models and Collective Phenomena in Membranes: Computer Simulation of Membrane Fusion, J. Polym. Sci. B: Polymer Physics, 41, 1441–1450, 2003. 42. Y. Kozlovsky and M. M. Kozlov, Stalk model of membrane fusion: solution of energy crisis, Biophys. J., 82, 882–895, 2002. 43. M. Müller and J. J. de Pablo, Computational approaches for the dynamics of structure formation in self-assembling polymeric materials, Annu. Rev. Mater. Sci., 43, submitted, 2013. 44. M. Müller and K. Ch. Daoulas, Speeding up intrinsically slow collective processes in particle simulations by concurrent coupling to a continuum description, Phys. Rev. Lett., 107, 227801, 2011. 18