VAGELIS HARMANDARIS International Conference on Applied Mathematics Heraklion, 16/09/2013

1. Simulations of Soft Matter under Equilibrium and Non-equilibrium Conditions VAGELIS HARMANDARIS International Conference on Applied Mathematics Heraklion, 16/09/2013

2. Outline • Introduction: Motivation, Length-Time Scales, Simulation Methods. • Multi-scale Particle Approaches: Atomistic and systematic coarse-grained simulations of polymers. • Application: Equilibrium polymeric systems. • Application: Non-equilibrium (flowing) polymer melts. • Conclusions – Open Questions.

3. COMPLEX SYSTEMS: TIME - LENGTH SCALES: • A wide spread of characteristic times: (15 – 20 orders of magnitude!) -- bond vibrations: ~ 10-15 sec -- dihedral rotations: 10-12sec -- segmental relaxation: 10-9 - 10-12 sec -- maximum relaxation time, τ1: ~ 1 sec (for Τ < Τm)‏

4. Hamiltonian (conserved quantity): Modeling of Complex Systems: Molecular Dynamics • Classical mechanics: solve classical equations of motion in phase space,Γ:=Γ(r, p). • In microcanonical (NVE) ensemble: The evolution of system from time t=0 to time t is given by : Liouville operator:

5. Modeling of Complex Systems: Molecular Dynamics • Various methods for dynamical simulations in different ensembles. • In canonical (NVT) ensemble: -- Langevin (stochastic) Thermostat -- Nose-Hoover thermostat: [Nosé 1984; Hoover, 1985]: add one more degree of freedom ζ.

6. -- stretching potential -- bending potential -- dihedral potential Van der Waals (LJ) Coulomb -- non-bonded potential Molecular Interaction Potential (Force Field) Molecular model: Information for the functions describing the molecular interactions between atoms. -- Potential parameters are obtained from more detailed simulations or fitting to experimental data.

7. MULTI-SCALE DYNAMIC MODELING OF COMPLEX SYSTEMS Atomistic MD Simulations: Quantitative predictions of the dynamics in soft matter. Limits of Atomistic MD Simulations (with usual computer power): -- Length scale: few Å – O(10 nm)‏ -- Time scale: few fs - O(1 μs) (10-15 – 10-6 sec) ~ 107 – 109 time steps -- Molecular Length scale (concerning the global dynamics): up to a few Me for “simple” polymers like PE, PB much below Me for more complicated polymers (like PS)‏ Need: - Simulations in larger length – time scales. - Application in molecular weights relevant to polymer processing. - Quantitative predictions. Proposed method: - Coarse-grained particle models obtained directly from the chemistry.

8. Systematic Coarse-Graining: Overall Procedure 1. Choice of the proper CG description. -- Microscopic (N particles) -- Mesoscopic (M “super particles”) -- Usually T is a linear operator (number of particles that correspond to a ‘super-particle’

9. Systematic Coarse-Graining: Overall Procedure 2. Perform microscopic (atomistic) simulations of short chains (oligomers) (in vacuum) for short times. 3. Develop the effective CG force field using the atomistic data-configurations. 4. CG simulations (MD or MC) with the new coarse-grained model. Re-introduction (back-mapping) of the atomistic detail if needed.

10. Effective (Mesoscopic) CG Interaction Potential (Force Field) CG Potential:In principle UCG is a function of all CG degrees of freedom in the system and of temperature (free energy): • CG Hamiltonian – Renormalization Group Map: • Remember: • Assumption 1:

11. r Bonded CG Interaction Potential Bonded Potential • Degrees of freedom:bond lengths (r), bond angles (θ), dihedral angles () Procedure: • From the microscopic simulations we calculate the distribution functions of the degrees of freedom in the mesoscopic representation,PCG(r,θ,). • PCG(r,θ, ) follow a Boltzmann distribution: • Assumption 2: • Finally:

12. q Non-bonded CG Interaction Potential: Reversible Work • Assumption 3: Pair-wise additivity Reversible work method [McCoy and Curro, Macromolecules, 31, 9362 (1998)] • By calculating the reversible work (potential of mean force) between the centers of mass of two isolated molecules as a function of distance: • Average < > over all degrees of freedom Γ that are integrated out (here orientational ) keeping the two center-of-masses fixed at distance r.

13. CG MD DEVELOPMENT OF CG MODELS DIRECTLY FROM THE CHEMISTRY APPLICATION: POLYSTYRENE (PS)‏ [Harmandaris, et al. Macromolecules, 39, 6708 (2006); Macromol. Chem. Phys. 208, 2109 (2007); Macromolecules 40, 7026 (2007); Fritz et al. Macromolecules 42, 7579 (2009)] 1) CHOICE OF THE PROPER COARSE-GRAINED MODEL • 2:1 model: Each chemical repeat unit replaced by two CG spherical beads (PS: 16 atoms or 8 “united atoms” replaced by 2 beads). • CG operator T: from “CHx” to “A” and “B” description. • Each CG bead corresponds to O(10) atoms. σΑ = 4.25 Å σB = 4.80 Å • Chain tacticity is described through the effective bonded potentials. • Relatively easy to re-introduce atomistic detail if needed. 2) ATOMISTIC SIMULATIONS OF ISOLATED PS RANDOM WALKS

14. CG MD Simulations: Structure in the Atomistic Level after Re-introducing the Atomistic Detail in CG Configurations. • Simulation data: atomistic configurations of polystyrene obtained by reinserting atomistic detail in the CG ones. • Wide-angle X-Ray diffraction measurements [Londono et al., J. Polym. Sc. B, 1996.] grem: total g(r) excluding correlations between first and second neighbors.

15. CG Free energy Atomistic Configuration CG Polymer Dynamics is Faster than the Real Dynamics PS, 1kDa, T=463K Free Energy Landscape • --CG effective interactions are softer than the real-atomistic ones due to lost degrees of freedom (lost forces). • This results into a smoother energy landscape. • CG MD: We do not include friction forces.

16. Time Scaling CG Polymer Dynamics – Quantitative Predictions CG dynamics is faster than the real dynamics. Time Mapping (semi-empirical) method: • Find the proper time in the CG description by moving the raw data in time. Choose a reference system. Scaling parameter, τx, corresponds to the ratio between the two friction coefficients. Time Mapping using the mean-square displacement of the chain center of mass • Check transferability of τx for different systems, conditions (ρ, T, P, …).

17. Polymer Melts through CG MD Simulations: Self Diffusion Coefficient • Correct raw CG diffusion data using a time mapping approach. • [V. Harmandaris and K. Kremer, Soft Matter, 5, 3920 (2009)] • Crossover regime: from Rouse to reptation dynamics. Include the chain end (free volume) effect. -- Rouse: D ~ M-1 -- Reptation: D ~ M-2 Crossover region: -- CG MD: Me ~ 28.000-33.000 gr/mol -- Exp.: Me ~ 30.0000-35.000 gr/mol -- Exp. Data: NMR [Sillescu et al. Makromol. Chem., 188, 2317 (1987)]

18. CG Simulations – Application: Non-Equilibrium Polymer Melts • Non-equilibrium molecular dynamics (NEMD): modeling of systems out of equilibrium - flowing conditions. • NEMD: Equations of motion (pSLLOD) simple shear flow • In canonical ensemble (Nose-Hoover) [C. Baig et al., J. Chem. Phys., 122, 11403, 2005] :

19. ux Primary y x x CG Simulations – Application: Non-Equilibrium Polymer Melts • NEMD: equations of motion are not enough: we need the proper periodic boundary conditions. • Steady shear flow: Lees-Edwards Boundary Conditions simple shear flow

20. CG Polymer Simulations: Non-Equilibrium Systems • CG NEMD - Remember: CG interaction potentials are calculated as potential of mean force (they include entropy). In principle UCG(x,T) should be obtained at each state point, at each flow field. Important question: How well polymer systems under non-equilibrium (flowing) conditions can be described by CG models developed at equilibrium? Method: [C. Baig and V. Harmandaris, Macromolecules, 43, 3156 (2010)] Use of existing equilibrium CG polystyrene (PS) model. • Direct comparison between atomistic and CG NEMD simulations for various flow fields. Strength of flow (Weissenberg number, Wi = 0.3 - 200) • Study short atactic PS melts (M=2kDa, 20 monomers) by both atomistic and CG NEMD simulations.

21. CG Non-Equilibrium Polymers: Conformations • Properties as a function of strength of flow (Weissenberg number) • Conformation tensor R • Atomistic cxx: asymptotic behavior at high Wi because of (a) finite chain extensibility, (b) chain rotation during shear flow. • CG cxx: allows for larger maximum chain extension at high Wi because of the softer interaction potentials.

22. CG Non-Equilibrium Polymers: Conformation Tensor • cyy, czz: excellent agreement between atomistic and CG configurations.

23. CG Non-Equilibrium Polymers: Dynamics • Is the time mapping factor similar for different flow fields? • [C. Baig and V. Harmandaris, Macromolecules, 43, 3156 (2010)] Translational motion • Purely convective contributions from the applied strain rate are excluded. • Very good qualitative agreement between atomistic and CG (raw) data at low and intermediate flow fields.

24. CG Non-Equilibrium Polymers: Dynamics Orientational motion • Rotational relaxation time: small variations at low strain rates, large decrease at high flow fields. • Good agreement between atomistic and CG at low and intermediate flow fields.

25. CG Non-Equilibrium Polymers: Dynamics • Time mapping parameter as a function of the strength of flow. • Strong flow fields: smaller time mapping parameter  effective CG bead friction decreases less than the atomistic one. Reason: CG chains become more deformed than the atomistic ones.

26. Conclusions • Hierarchical systematic CG models, developed from isolated atomistic chains, correctly predict polymer structure and dimensions. • Time mapping using dynamical information from atomistic description allow for quantitative dynamical predictions from the CG simulations, for many cases. • Overall speed up of the CG MD simulations, compared with the atomistic MD, is ~ 3-5 orders of magnitude. • System at non-equilibrium conditions can be accurately studied by CG NEMD simulations at low and medium flow fields. • Deviations between atomistic and CG NEMD data at high flow fields due to softer CG interaction potentials.

27. Challenges – Current Work • Estimation of CG interaction potential (free energies): Check – improve all assumptions • Ongoing work with M. Katsoulakis, D. Tsagarogiannis, A. Tsourtis • Quantitative prediction of dynamics based on statistical mechanics • e.g. Mori-Zwanzig formalism (Talk by Rafael Delgado-Buscalioni) • Parameterizing CG models under non-equilibrium conditions • e.g. Information-theoretic tools (Talk by Petr Plechac) • Application of the whole procedure in more complex systems • e.g. Multi-component biomolecular systems, hybrid polymer based nanocomposites • Ongoing work with A. Rissanou

28. ACKNOWLEDGMENTS Prof. C. Baig [School of Nano-Bioscience and Chemical Engineering, UNIST University, Korea] Funding: ACMAC UOC [Regional Potential Grant FP7] DFG [SPP 1369 “Interphases and Interfaces ”, Germany] Graphene Research Center, FORTH [Greece]

29. EXTRA SLIDES

30. APPLICATION: PRIMITIVE PATHS OF LONG POLYSTYRENE MELTS • Describe the systems in the levels of primitive paths [V. Harmandaris and K. Kremer, Macromolecules, 42, 791, (2009)] • Entanglement Analysis using the Primitive Path Analysis (PPA) method [Evereaers et al., Science 2004, 303, 823]. PP PS configuration (50kDa) CG PS configuration (50kDa) • Calculate directly PP contour length Lpp,, tube diameter: -- PP CG PS: Ne ~ 180 ± 20 monomers

31. CALCULATION OF Me in PS: Comparison Between Different Methods • Several methods to calculate Me: broad spread of different estimates [V. Harmandaris and K. Kremer, Macromolecules, 42, 791 (2009)]

32. MESOSCOPIC BOND ANGLE POTENTIAL OF PS Distribution function PCG(θ,T)‏ CG Bending potential UCG(θ,T)‏

33. CG Simulations – Applications: Equilibrium Polymer Melts • Systems Studied: Atactic PS melts with molecular weight from 1kDa (10 monomers) up to 50kDa (1kDa = 1000 gr/mol). • NVT Ensemble. • Langevin thermostat (T=463K). • Periodic boundary conditions.

34. STATIC PROPERTIES : Radius of Gyration RG

35. SMOOTHENING OF THE ENERGY LANDSCAPE Qualitative prediction: due to lost degrees of freedom (lost forces) in the local level  Local friction coefficient in CG mesoscopic description is smaller than in the microscopic-atomistic one Rouse: Reptation: CG diffusion coefficient is larger than the atomistic one

36. Time Mapping Parameter: Translational vs Orientational Dynamics