Simulation studies of soft matter: generic statistical properties and chemical details

The relation between atomistic structure, architecture, molecular weight and material properties is a basic concern of modern soft material science. This by now goes far beyond standard properties of bulk materials. A typical additional focus is on surface or interface aspects or on the relation between structure and function in nanoscopic molecular assemblies. This all implies a thorough understanding on many length and correspondingly time scales ranging from (sub)-atomic to macroscopic. At this point computer simulations are playing an increasingly important, if not the central role. Traditionally simulations have been separated in two main groups, namely simplified models to deal with generic or universal aspects of polymers, i.e. critical exponents, and those employing classical force field simulations with (almost) all atomistic detail, i.e. for the diffusion of small additives in a small “sample”. Still characteristic problems, which require huge systems and/or long times in combination with a chemistry specific model, cannot be tackled by these methods alone. More recently with the development of scale bridging or multi scale simulation techniques, these different approaches have been combined into an emerging rather powerful tool. It is the purpose of this contribution to give a few examples of how such an approach can be used to understand specific material properties.


Introduction
Synthetic and biological macromolecules, (bio-)membranes, colloidal suspensions etc. are classified under the general term soft matter. Soft means that the characteristic energy densities, which to a first approximation are close to the elastic constants, are several orders of magnitude smaller than for conventional "hard" matter. The relevant energy scale is the thermal energy k B T . This goes in hand with the fact that the conformational entropy of a macromolecule is of the same order as the intermolecular energy and their interplay usually determines the properties.
Macromolecules/polymers have been extensively studied on the basis of rather simple generic models. These ideas date back to the seminal work of Flory [1,2] and others and the subsequent solid theoretical basis and link to critical phenomena by de Gennes [3], who showed that all polymers follow the same universal scaling laws and that one can understand qualitatively and semi quantitatively their behavior on the basis of scaling theories. I.e. the mean squared extension of a polymer, R 2 (N ) , generally follows the scaling law R 2 (N ) = AN 2ν , where ν is a universal critical exponent. In so called good solvent a e-mail: kremer@mpip-mainz.mpg.de ν 0.6(d = 3), while it is 1/2 in a polymer melt. Similar relations hold for many other properties, including dynamics [3,4]. Though the power laws are universal, scaling theory says nothing about the prefactor A, which is determined by the chemical structure of the polymer and the interaction with its surrounding and can vary significantly. When it comes to mixtures, molecular assemblies, functional systems etc. one is often not fully reaching the scaling regime and the combination of generic and chemistry specific information is absolutely conditional for any theoretical prediction. Figure 1 illustrates the different scales and the different characteristic levels of simulation models. Many phenomena in biology, chemistry, and materials science involve processes occurring on atomistic length and time scales, which affect structural and dynamical properties on scales far beyond atomistic ones. Because it is infeasible (and most often undesirable!) to run computer simulations of very large systems with atomistically detailed models, mesoscale (coarse grained, cg) models are being developed, which stay reasonably close to the chemical structure. Reintroduction of chemical details allows to study atomistically detailed processes in various windows of the cg trajectory [5][6][7][8][9][10]. Such a procedure requires the (high-and low-) resolution models to be (pairwise) structurally consistent. However, since many high resolution  states correspond to one low resolution cg-state the conformational statistics on both levels of resolution should be the same if analyzed on the cg level.
In this paper we will emphasize a structure-based coarse graining whose aim is to allow for structurebased scale hopping [11]. Alternative coarse graining approaches [12][13][14] as well as approaches, which go much further and map the whole chain to one ellipsoidal particle [15,16] or just a soft sphere [17] will not be discussed here. For a general overview of advanced computer simulation techniques I refer to references [18,19].

Structure-based coarse graining approach
The starting point is a mapping scheme, which simplifies the interactions significantly, while staying close enough to the chemical structure to allow for switching between scales [6,7,10,11]. Figure 2 gives an example for polycabonate (BPA-PC) and two recently tested mapping schemes for polystyrene (PS). To generate the interaction potentials for the cg system we make the ansatz to strictly separate the parametrization of bonded and nonbonded interactions, based on the assumption that the total potential energy can be separated into a bonded/covalent, U cg B , and nonbonded, U cg NB . For U cg B of the coarse grained (cg) bond and torsion angles and bond lengths we sample distribution functions of the underlying atomistic model in vacuum by a stochastic simulation, i.e. Monte Carlo. To avoid any double counting of interactions with respect to the later treated nonbonded interactions, correlations along the backbone are only considered up to the distance needed to generate the appropriate distributions. The underlying atomistic potentials preferably originate from quantum calculations, thus are not linked to any parametrization based on experiments, which typically include already further interactions. The distributions in general are characterized by the cg bond lengths r, bond angles θ, and torsions φ, i.e., P (r, θ, φ, T ). In many cases several different bond and torsion angles have to be considered. For simplicity, one usually assumes that the distributions decouple, namely P (r, θ, φ, T ) = P (r, T ) P (θ, T )P (φ, T ). This assumption generally does not hold and the introduced error strongly depends on details of the choice of the cg model [20], which of course has to be checked carefully. Boltzmann inversion of the distributions yield temperature dependent bonded cg interaction potentials, i.e. for the bond length The temperature dependence concerns not only the prefactor k B T but also the distributions P themselves. Strictly speaking U cg can only be applied at the temperature they were derived at. For apolar, amorphous polymers nonbonded interactions are usually sufficiently described by an excluded volume potential, where the parameters are either estimated by group contribution methods or from the van der Waals volumes of the beads. For BPA-PC the latter was sufficient, leading to simple repulsive Lennard-Jones potentials (WCA potential) [22]. More systematically radial distributions functions of small all-atom and cg chains under appropriate melt or solution conditions are used to derive an appropriate tabulated potential by e.g. the iterative Boltzmann inversion method [13,21]. For BPA-PC a characteristic ratio R 2 G /N ≈ 37Å 2 , where R G is the radius of gyration and N the number of repeat units, is found and compares very well with neutron scattering experiments [9] and allowed us to study structural as well as dynamic properties of BPA-PC melts for systems up to 200 chains up to 20 entanglement molecular weights (N e ≈ 5-6, N = 120), corresponding to system sizes of (100 nm) 3 [23,24]. In a similar way PS can be studied and a similar level of agreement is reached [20] giving R 2 G /N ≈ 16Å 2 at T = 463 K.

Inverse mapping
In the previous section we gave a very short account on the structure based coarse graining procedure for amorphous polymers. Since the cg models are still fairly close to the original chemical structure we can also make the step back and reintroduce the chemical details of the chains. That then allows for a detailed comparison to microscopic experiments. It also opens the way to the anticipated scale hopping or adaptive resolution simulation. There are several options to reintroduce chemical details [6,7,11,22,24,25]. If the polymers consist of reasonably rigid (all-atom) fragments it is sufficient to fit these rigid all-atom units onto the corresponding cg chain segment coordinates. The atomistic fragments are taken from a pool of structures, that correctly reflect the statistical weight of those degrees of freedom (DOFs), i.e., ring flips, that are not resolved in the cg description. Then a very short equilibration run of a few ps is performed, where the individual atoms typically move a distance less than 1.5Å. This is to be compared to the typical strand-strand distance of about 5Å in a polymer melt or to the size of a monomer, which in the case of BPA-PC is about 11Å. Thus this is just a very local rattling. The so derived structures compare very well to scattering experiments as shown in Figure 3. Here we compare BPA-PC data of all atom melts, derived from a cg simulation to neutron scattering experiments and PS united atom melts (all H atoms are included in "carbon super atoms") also derived from cg simulations to wide angle X-ray scattering.
Together with the characteristic ratio this indicates that the obtained all-atom melts truly represent the experimental system, supporting the overall consistency of the approach. Note that this procedure does not require any fitting to experiment for the bonded part and, for the present case, not even for the nonbonded part! In a first application to phenol diffusion in BPA-PC [26] the polymer matrix was generated by the above procedure to host the added small molecules.

Dynamics
The structure based mapping scheme fixes the length scaling by its very construction. This is not the case for dynamics. cg simulations are usually analyzed in Lennard Jones units , σ and τ . Using realistic values for these units (mass of a cg unit in the model etc.) one can determine the time scale of the cg simulation. Giving the same mass to all cg beads results in 1τ = 1.7 ps for BPA-PC at the process temperature of 570 K and 1τ = 1.0 ps for PS at 463 K respectively. While these are the intrinsic time scales of the cg models, it is not the correct time scale for the underlying chemical system [23,24]. Although on the cg level scaling properties of the dynamics, such as Rouse behavior for short chain melts and reptation behavior of long chain melts are well reproduced and the entanglement length is predicted very precisely, the time scale has yet to be determined. The reason for that originates from the fact, that the polymer melt dynamics is dominated by the chain chain friction, which results from local interactions as well as packing and topological constraints. Because of that it is not possible to calculate the time scaling from the potential functions of the simulations. A pragmatic direct way is to match mean square displacements of short chain cg melt simulations to those of all or united atom chain melts. Doing this not only N5, 290 K N5, 570 K N20, 290 K N20, 570 K experiment Fig. 3. Comparison of experimentally obtained structure information to results from the simulations after reintroducing the details. (a) Structure factor of a polycarbonate melt compared to neutron scattering (upper panel, after Ref. [24]). (b) Inter chain radial distribution functions for PS compared to X-ray scattering results (lower panel, after Ref. [27]).
gives the time scaling but also shows, down to what scale the cg description follows the appropriate molecular motion. In the case of BPA-PC at 570 K we find 1τ ≈ 30 ps and that the cg displacements properly represent atomistic trajectories down to distances of about 10Å, roughly the size of a monomer. Combining all this means that in our case we now have systems with about 800 000 atoms for times of up to 5 × 10 −5 s at hand. Along these lines recently we also investigated a variety of problems related to polystyrene [20,27] and find 1τ ≈ 7.7 ps. Taking all this gives us confidence that the structure based coarse graining approach yields solid information on the statics and dynamics of complicated polymeric melts. Estimates of the speedup for the current simulations are i.e. around 10 4 for BPA-PC compared to brute force all atom simulations.

528
The European Physical Journal B

Extensions and outlook
This ansatz has been used to study a number of different physical problems, which otherwise would have been too time demanding. The extensive study of entangled polycarbonate melts [23] has already been mentioned before. Already somewhat earlier we used the inversely mapped cg structures of polycarbonate as hosts for a few phenol molecules, in order to study their diffusion. There we found that phenol represents an interesting case of activated diffusion in between the small atom or small molecule hopping motion and the more matrix fluctuation mediated diffusion of large molecules, such as added oligomers [26]. Also the observed vacancy structure compared well to positron annihilation experiments. By combining united atom simulations with cg simulations and NMR experiments we recently studied the diffusion of ethyl benzene in a PS matrix at a variety of temperatures and concentrations. Beyond that the interaction of polymers with specific surfaces is of special interest. While dispersion forces lead to an overall unspecific attraction of the polymers to the surface, the resulting changes of the melt conformations are well described by rather idealized theories. In the case of a metal surface, such as Ni, however, specific interactions of fragments of the chains might lead to rather different morphologies, as shown for the example of BPA-PC [28]. There the above cg scheme was combined with an ab initio density functional MD study of specific surface interactions. For this the monomer had to be cut into pieces and the interaction of the submolecules benzene and phenol, carbonic acid and propane with a Ni(1, 1, 1) surface was studied [29]. Because both carbonic acid and propane are strongly repelled by the Ni surface the strong attraction experienced by benzene is screened within the chain. This leads to the situation where only the terminal benzenes can adsorb. The comparison of different chain ends shows that the melt morphology close to the surface for phenolic chain ends is significantly different from the generic dispersion force dominated structure [29]. This has also been extended to study the sliding of BPA-PC melts under shear with and without additives [30][31][32] and at a surface defect [33].
The above described examples all treat a simulation on one specific level of detail and then compare them. What instead would be highly desirable is an approach, where (two) levels of resolution can be applied at the same time within one simulation in different areas and where both regimes are in true equilibrium with each other, allowing for a free exchange of molecules. Then one could run a simulation on a lower level of detail and zoom in on-the-fly, whenever "something interesting" happens. Recently, we proposed the particle-based Adaptive Resolution Scheme (AdResS) [34,35] for hybrid all-atom/coarsegrained molecular dynamics (MD) simulations that fulfils the above requirement. In a first application a model of liquid water [36] was studied by that approach. These however are only very first steps and a detailed description would be beyond the scope of this short contribution.  Over the years I have enjoyed the collaboration with many colleagues and students on the topics described in this short overview. Out of them I would like to especially thank C.F.