1 / 101

VAGELIS HARMANDARIS

Hierarchical Multi-scale Modeling of Polymers under Equilibrium and Non-equilibrium conditions: Computational and Mathematical Aspects. VAGELIS HARMANDARIS ACMAC Workshop: Coarse-graining of many-body systems: analysis, computations and applications Crete, 01/07/2011. Outline.

edie
Download Presentation

VAGELIS HARMANDARIS

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Hierarchical Multi-scale Modeling of Polymers under Equilibrium and Non-equilibrium conditions: Computational and Mathematical Aspects VAGELIS HARMANDARIS ACMAC Workshop: Coarse-graining of many-body systems: analysis, computations and applications Crete, 01/07/2011

  2. Outline • Introduction: Motivation, Length-Time Scales. • Multi-scale particle approaches: Atomistic and systematic coarse-grained simulations of polymers. • Application I: Equilibrium polymers. • Application II: Non-equilibrium (flowing) polymers. • Application III: Polymer/solid interfaces. • Conclusions – Open questions.

  3. -- Clever Materials • Molecular Electronics, Carbon Nanotubes COMPLEX SYSTEMS -- Broad spectrum of systems, applications, length-time scales. • Systems -- polymers -- biological macromolecules (cell membrane, DNA, lipids)‏ -- colloids -- hybrid interfacial systems -- ... • Applications --- Nanotechnology (materials in nano-dimensions), biotechnology (drug release, … etc)‏

  4. Complex Macromolecular 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 (σε Τ < Τm)‏

  5. Β) description at microscopic (atomistic) level C) description at mesoscopic (coarse-grained) level D) description at continuum level Α) description at quantumlevel Hierarchical Modeling of Molecular Materials • Due to wide range of characteristic lengths - times, several simulation methods that describe length and time scales have been developed:

  6. Atomistic Molecular Dynamics Method • Molecular Dynamics (MD)[Alder and Wainwright, J. Chem. Phys., 27, 1208 (1957)] • Assume a microscopic (atomistic) system of N interacting particles. Each particle has position ri = (xi,yi, zi). In 3D there are 6N degrees of freedom (3N positions + 3N velocities). Phase space is R6N. • Particles interact with a very complicated and in principle unknown function of 3N variables: U(r1,r2,..., rN). • The time evolution of each particle is (Newton’s eqs. of motion): • System of 6N 1st order PDE’s. N ~ 104 – 105 particles

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

  8. Polymer Dynamics through Atomistic MD Simulations • Atomistic MD simulations: capable of predicting accurately polymer dynamics; at least for simple polymers! Systems: various polymer lengths from C40, to C300 (M=4800) for times up to 200 ns. -- T = 450 K > Tm -- Equilibration: using a powerful MC algorithm (Theodorou et al, 1999) -- TraPPE model (Siepmann et al, 1998) -- Multiple Time Step algorithm (Martyna et al., 1996) Polyethylene, T=450K [Harmandaris et al. Macromolecules, 1998; 2000; 36, 1376, 2003]

  9. MULTI-SCALE MODELING OF COMPLEX SYSTEMS Atomistic MD Simulations: Allows for quantitative predictions (real numbers!) 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 Ne for “simple” polymers like PE, PB much below Ne for more complicated polymers (like PS)‏ Need: - Simulations in larger length – time scales. - Application in molecular weights relevant to polymer processing. One solution: - Coarse-grained particle models obtained directly from the chemistry.

  10. CG METHODS: MODELS DEVELOPED DIRECTLY FROM THE CHEMISTRY [Luybarchev, 1997; Tschöp et al, 1998; F.M. Plathe, 2001; Voth, 2002; Harmandaris et al, 2006; Peter et al. 2008] • Goal: predict properties directly from first-principles. • Integrate out some degrees of freedom as one moves from finer to coarser scales. • Main Idea: Information from the atomistic level is used in the development of the CG force field. Microscopic (atomistic) simulations CG Force field Mesoscopic CG simulations • Directly linked to the chemistry. • Allow for quantitative predictions of material properties.

  11. Systematic Coarse-Graining: Overall Procedure to Develop CG Models Directly from the Chemistry 1. Choice of the proper CG description. -- number of atoms that correspond to a ‘super-atom’ (coarse grained bead)‏ 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.

  12. DEVELOPMENT OF CG MODELS DIRECTLY FROM THE CHEMISTRY APPLICATION: POLYSTYRENE (PS)‏ [VH, et al. Macromolecules, 39, 6708 (2006); Macromol. Chem. Phys. 208, 2109 (2007); Macromolecules 40, 7026 (2007); Fritz et al. Macromolecules 42, 7579 (2009)] Choice of the Proper CG Variables: Depends on the Physical Problem • 2:1 model: Each chemical repeat unit replaced by two 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.

  13. 2) ATOMISTIC SIMULATIONS OF PS OLIGOMERS • TraPPE UA model[Siepmann and coworkers, 2000]. • All-atom model [Plathe, 1996]. A1) Single Polymer Chain A2) Potential of Mean Force between two Oligomers • All possible PS tacicities: isotactic (RS/SR), syndiotactic (RR/SS). • Used to develop the effective CG force field. B) Long Atomistic Simulations of Short Chain Polymer Melts: • very long (up to 1.0 μs) atomistic MD runs for very short PS melts (up to 3kDa) are performed. • Used to compare with CG data.

  14. CORSE-GRAINED MOLECULAR DYNAMICS • CG representation: • Instead of N atomistic particles we have M CG beads (with M < N). • Each particle has coordinates qi = (xi,yi, zi) . There are 3M CG degrees of freedom, instead of 3N atomistic ones. • Particles interact with a general and in principle unknown function UCG(q1,q2,..., qM)= UCG(Q). • Evolution equation (Brownian motion) of the CG-mesoscopic system. • Important (still open) mathematical problem: • How do we accurately derive UCG(Q) from UAT(r)? • CG Hamiltonian – Renormalization Group Map:

  15. DERIVATION OF THE EFFECTIVE MESOSCOPIC FORCE FIELD USING ATOMISTIC DATA: CG Potential:In principle UCG(q) is a function of all degrees of freedom in the system. Assumption: Decomposition of UCG(q) in bonded and non-bonded part: • BONDED POTENTIAL • Degrees of freedom: CG bond lengths (rCG), bond angles (θCG), dihedral angles (φCG)‏. r • Question: How do we calculate UbondedCG (rCG,θCG,CG , …)?

  16. IDEA: CG Particles follow Boltzmann distributions Procedure: From the microscopic simulations we calculate the distribution functions of the degrees of freedom in the CG representation, PCG(rCG,θCG,φCG). Then by inverting the Boltzmann relation we calculate effective potentials (free energies) in the CG representation, UbondedCG (rCG,θCG,CG , …). • Assumption: • Finally:

  17. NONBONDED INTERACTION PARAMETERS: REVERSIBLE WORK • CG Hamiltonian – Renormalization Group Map: • 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.

  18. q NONBONDED INTERACTION PARAMETERS: REVERSIBLE WORK • Calculate “reversible work” using a numerical method (eg. MC or MD). • Main assumptions: • (A) Neglect many-body effect. Exact for the gas phase. • (B) Chain effect is not described. In our case CG particles belong to a macromolecule. Solution: Use of conditional reversible work [Fritz et al, 2009; E. Brini et al. 2011] • Main idea: Use instead two very short chains and keep constant the distance between the center-of-mass of only the two target CG (e.g. A, B) particles. • Calculate the PMF including all atomistic interactions, • Calculate the PMF with all atomistic interactions excluding the A-B ones. • Effective CG interaction is:

  19. CG NONBONDED EFFECTIVE INTERACTION POTENTIAL • CG effective potential calculated by the reversible work method using short chains.

  20. 4) CG MOLECULAR DYNAMICS SIMULATIONS • System: Atactic PS melts with molecular length from 1kDa (10 monomers) up to 10kDa (1kDa = 1000 gr/mol, T=463K). • Method: MD with a Langevin thermostat • Friction Force: with friction coefficient Γ = 1.0/τ • Random force: (fluctuation-dissipation theorem)

  21. 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.

  22. BACK-MAPPING PROCEDURE: FROM MESOSCOPIC DESCRIPTION TO ATOMISTIC DETAIL -- In many cases we need microscopic configuration. • Goal: from the CG configuration obtain a microscopic one. • Note: One CG configuration corresponds to many microscopic ones No unique solution !!

  23. Coarse - Grained 2:1 PS model A B BACK-MAPPING PROCEDURE: FROM MESOSCOPIC DESCRIPTION TO ATOMISTIC DETAIL [Harmandaris, et al. Macromolecules, 39, 6708 (2006)] • Example: PS CG melt -- Initial configuration: coarse-grained PS 2:1 melt -- Final configuration: microscopic (atomistic) PS melt

  24. BACK-MAPPING PROCEDURE: FROM MESOSCOPIC DESCRIPTION TO ATOMISTIC DETAIL [Harmandaris, et al. Macromolecules, 39, 6708 (2006)] • Back-mapping in 3 steps: • STEP 1: Geometrical Constructionof Each Atomistic Chain Using the CG Configurations • Every monomer is constructed with a numerical method using the CG points of the current monomer and of the neighboring monomers. • Minimization of an objective function F: -- Qij: equilibrium distance between atom i and center-of-mass of CG bead j. -- Quasi-Newton algorithm.

  25. BACK-MAPPING PROCEDURE: FROM MESOSCOPIC DESCRIPTION TO ATOMISTIC DETAIL Total Energy Minimization Scheme • STEP 2: Very Short (~20 ps) MD Run to Obtain Realistic Configurations • STEP 3:

  26. STATIC PROPERTIES :LOCAL STRUCTURE • radial distribution function gn(r): describe how the density of surrounding matter varies as a distance from a reference point. • pair radial distribution function g(r)=g2(r): gives the joint probability to find 2 particles at distance r. Easy to be calculated in experiments (like X-ray diffraction) and simulations. • choose a reference atom and look for its neighbors:

  27. VALIDATION OF THE BACK-MAPPING PROCEDURE • Compare structure in atomistic configurations (of a 1kDa system) obtained from: • Very long (0.5 μs) “pure” atomistic MD runs. • Short-fast CG MD runs: reinsertion of the atomistic detail in CG configurations. Intramolecular g(r) Intermolecular CH2 – CH2 g(r)‏

  28. STRUCTURE IN THE ATOMISTIC LEVEL USING THE BACK-MAPPING PROCEDURE • 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.

  29. CG Free energy Atomistic Configuration CG Dynamics is Faster than the Microscopic (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.

  30. 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

  31. Friction – Noise Forces • --CG effective interactions are softer than the real-atomistic ones due to lost degrees of freedom (lost forces). • NOTE: In CG MD, we do not (properly) include friction forces. Generalized Langevin equation: We use a simple Langevin equation where friction force is acting only as a thermostat:

  32. Time Scaling CG Polymer Dynamics – Quantitative Predictions • Semi-empirical method: Find the proper time in the CG description by moving the raw data in time. Scaling parameter, corresponds to the ratio between the two friction coefficients (scalar parameter in homogeneous systems). • Perform long microscopic (atomistic) simulations in a reference system. Time mapping at that specific state point: calculate τx for this system! Time Mapping using the mean-square displacement of the chain center of mass • Check transferability of τx for different systems, conditions (ρ, T, P, …).

  33. Time Mapping Parameter: Translational vs Orientational Dynamics

  34. CG Polymer Dynamics: 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)]

  35. Polymer Segmental Dynamics: Effect of temperature and Pressure • Study segmental dynamics as a function of temperature and pressure. • [V. Harmandaris, et al., Macromolecules 44, 393 (2011)] • Calculate auto-correlation time functions of bond order parameters:

  36. Polymer Segmental Dynamics: Effect of temperature and Pressure • Get segmental relaxation times from MD simulations: Compare directly with dielectric experiments.

  37. CG Simulations – Application II: Non-Equilibrium Polymer Melts • Non-equilibrium molecular dynamics (NEMD): modeling of systems out of equilibrium - flowing conditions. • NEMD: Equations of motion in canonical (NVT) ensemble [C. Baig et al., J. Chem. Phys., 122, 11403, 2005] simple shear flow Lees-Edwards Boundary Conditions

  38. ux Primary y x x CG POLYMER DYNAMICS: NON-EQUILIBRIUM SYSTEMS • NEMD: equations of motion are not enough: we need the proper periodic boundary conditions. • Steady shear flow: Lees-Edwards Boundary Conditions

  39. 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? Idea: Use an existing equilibrium CG polymer model under non-equilibrium conditions. • Direct comparison between atomistic and CG NEMD simulations for various flow fields. Strength of flow (Weissenberg number, Wi = 0.3 - 200) • Short atactic PS melts (M=2kDa, 20 monomers) are studied by both atomistic and CG NEMD simulations.

  40. CG Non-Equilibrium Polymers: Conformations • Properties as a function of strength of flow (Weissenberg number) • Conformation tensor: -- R is the end-to-end vector • 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.

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

  42. CG Non-Equilibrium Polymers: Dynamics • Is the time mapping factor similar for different dynamical quantities? 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.

  43. 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.

  44. 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.

  45. Application III: Polymer Nanocomposites • Hybrid Solid/Liquid interfacial systems. • Molecule/solid interface is of particular importance. Need an accurate molecule/surface interaction potential. Example: Polymer/Graphite System.

  46. Polymer – Graphite Interaction: Analytic Estimation Graphite substrate [Steele’s method, 1973]: • Takes into account the crystallographic structure of the substrate. -- Assume perfect crystal of graphite periodic in x,y andz<0 directions. -- Each polymer atom interacts with every graphite atom with a via a Lennard Jones potential: σps= 3.304Å, εps= 0.06839kcal/mol -- The total interaction energy for a polymer atom is the summation of every graphite atom/polymer atom interaction.

  47. Atomistic Molecular Dynamics Simulations: Applications Polyethylene/Graphite Interfacial Systems:Model of Polymer Nanocomposites [Harmandaris, et al. Macromolecules, 38, 5780 (2005); 38, 5796 (2005)] Model System: A thin polymer film adsorbed onto a solid graphite surface on one side and exposed to vacuum on the other • Various polymer lengths from C40, to C400 (M=4800) for times up to 100 ns. • T = 450 K > Tm • Equilibration: using a powerful MC algorithm (Theodorou et al, 1999) • TraPPE model (Siepmann et al, 1998) • Multiple Time Step (Martyna et al., 1996) z

  48. INTERFACIAL SYSTEMS: -- Transport diffusion coefficient at different distances from the solid surface: In overall: Effect of surface on density ~ 18 Å, conformations ~ 2 RG conformational dynamics ~ 6 Å global dynamics ~ 6-7 RG

  49. Polymer/Metal Interfaces • Molecule/Metal interaction cannot be estimated analytically! • Use detailed DFT calculations of a single molecule adsorbed onto the metal surface. Test example: Benzene/Au systems [K. Johnston and VH, J. Phys. Chem. C to be published] • Obtain: Molecule/Metal (Adsorption) energies of a single molecule for various: • (A) molecule conformations and • (B) distances from the surface.

  50. Density Functional Theory Density Functional Theory:

More Related