Studies of Nano, Chemical, and Biological Materials by Molecular Simulations

1. Studies of Nano, Chemical, and Biological Materials by Molecular Simulations Yanting Wang Institute of Theoretical Physics, Chinese Academy of Sciences Institute of Theoretical Physics, Chinese Academy of Sciences Beijing, China September 25, 2008

2. Atomistic Molecular Dynamics Simulation • Solving Newton’s Equations of Motion. • Empirical force fields are determined by fitting experimental results or data from first principles calculations • Quality of empirical force fields has big influence on simulation results • Capable of simulating up to millions of atoms (parallel computing)

3. Quantifying Condensed Matter Structures • Bond-Orientational Order Parameters • Capture the symmetry of spatial orientation of chemical bonds • Non-zero values for crystal structures • 0 for liquid • Radial Distribution Function g(r) • Appearance probability of other atoms with respect to a given atom • Discrete values for solids • Continuous waves for liquids • 1 for ideal gas (isotropic structure)

4. Some Applications of Gold Nanomaterials Electronic lithography Ion detection Molecular electronics Larger Au particles change color Gold nanowires Chemical etching J. Zheng et al., Langmuir 16, 9673 (2000) R. F. Service, Science 294, 2442 (2001) S. O. Obare et al., Langmuir 18, 10407 (2002) • Both size and shape are important in experiments!

5. Thermal Stability of Low Index Gold Surfaces • Stable gold interior: FCC structure • Thermal stability of surface: {110} < {100} < {111}

6. Stability of Icosahedral Gold Nanoclusters* Simulated annealing from a liquid • Empirical glue potential model • Constant T molecular dynamics (MD) • From 1500K to 200K with T=100K, and keep T constant for 21 ns • thousands of atoms Cooling • Strained FCC interior • All covered by stable {111} facets Liquid at T=1500K Icosahedron at T=200K • Mackay Icosahedron with a missing central atom • Asymmetric facet sizes * Y. Wang, S. Teitel, C. Dellago Chem. Phys. Lett. 394, 257 (2004) * Y. Wang, S. Teitel, C. Dellago J. Chem. Phys. 122, 214722 (2005)

7. First-Order Like Melting Transition Heat to melt • Keep T constant for 43 ns • T = 1075K for N = 2624 • Magic and non-magic numbers • Cone algorithm* to group atoms into layers Potential energy vs. T Surface Interior Sub-layers • First-order like melting transition * Y. Wang, S. Teitel, C. Dellago J. Chem. Phys. 122, 214722 (2005)

8. No Surface Premelt for Gold Icosahedral Nanoclusters Interior Surface N = 2624 • Interior keeps ordered up to melting temperature Tm • Surface softens but does not melt below Tm

9. Surface Atoms Diffuse Below Melting Mean squared displacements (average diffusion) N = 2624 • All surface atoms diffuse just below melting • Surface premelting?

10. “Premelt” of Vertices and Edges but not Facets Movement t=1.075ns Movement 4t Average shape Mechanism • Vertex and edge atoms diffuse increasingly with T • Facets shrink but do not vanish below Tm=1075 K • Facet atoms also diffuse below Tm because the facets are very small !

11. Conclusions • First-order like melting transition for gold nanoclusters with thousands of atoms • Very stable {111} facets result in good thermal stability of icosahedral gold nanoclusters • Vertex and edge “premelt” softens the surface but no overall surface premelting

12. Very Small Gold Nanoclusters? • Smaller gold nanocluster has more active catalytic ability • Debate if very small gold nanoclusters (< 2 nm ) are solid or liquid • 54 gold atoms (only two layers) • Not an icosahedron All surface atoms are on vertex or edge!

13. Smeared Melting Transition for N = 54* • Heat up sequentially • timestep 2.86 fs • 108 steps at each T • Easy to disorder due to less binding energy • Melting transition from Ts≈ 300 K to Te ≈ 1200 K Te Ts Heat capacity Average potential energy per atom * Y. Wang, S. Rashkeev J. Phys. Chem. C 113, 10517 (2009).

14. Snapshots at Different Temperatures • Both layers premelt below 560 K • No inter-layer diffusion below 560 K

15. Ti Ti Td Inter- and Intra- Layer Diffusion Moved atoms: moving to the other layer at least once at each temperature Atomic self diffusion starts at Td≈ 340 K Inter-layer diffusion starts at Ti≈ 560 K Liquid crystal-like structure between 340 K and 560 K

16. More Layers in Between: Approaching First-Order Melting Transition* Onset Temperature Ts and Complete Temperature Te of Melting Transition, Self Diffusion Temperature Td, and Interlayer Diffusion Temperature Ti • Melting temperature region narrows down for more layers • Only two-layer cluster has intra-layer diffusion first * Y. Wang, S. Rashkeev J. Phys. Chem. C 113, 10517 (2009).

17. Conclusions • Smeared melting transition for two-layer gold nanocluster • Mechanism consistent with icosahedral gold nanoclusters • Liquid-crystal like partially melted state for two-layer gold nanocluster: intra-layer diffusion without inter-layer diffusion • Approaching well-defined first-order melting transition for gold nanoclusters with more layers • Very small gold nanoclusters have abundant phase behavior that can not be predicted by simply extrapolating the behavior of larger gold nanoclusters

18. T=5K T=515K T=1064K T=1468K Thermal Stability of Gold Nanorods* Experimental model • Increasing total E continuously to mimic laser heating Pure FCC interior Z. L. Wang et al., Surf. Sci. 440, L809 (1999) Two steps * Y. Wang, C. Dellago J. Phys. Chem. B 107, 9214 (2003).

19. Surface-Driven Bulk Reorganization of Gold Nanorods* Surface Second sub layer Cross sections • Temperature by temperature step heating • Minimizing total surface area • Surface changes to all {111} facets • Interior changes fcc→hcp→fcc by sliding planes, induced by surface change • Interior fcc reorients Yellow: {111} Green: {100} Red: {110} Gray: other Yellow: fcc Green: hcp Gray: other * Y. Wang, S. Teitel, C. Dellago Nano Lett. 5, 2174 (2005).

20. Conclusions • Thermal stability of gold nanoclusters and gold nanorods is closely related to specific surface structures (not only surface stress matters) • Shape change of gold nanorods comes from the balance between surface and internal free energetics

21. Assuming central pairwise effective forces • Minimizing force residual Benifit: maller numbers of degrees of freedom and faster dynamics • Well rebuild structural properties • Can eliminate some atoms at CG level • Does NOT consider transferability! Multiscale Coarse-Graining (MS-CG) Method* to Rigorously Build CG Force Fields from All-Atom Force Fields • Pioneer work by Dr. Sergey Izvekov with block-averaging • Theory by Prof. Will Noid (Penn State U), Prof. Jhih-Wei Chu (UC-Berkeley), Dr. Vinod Krishna, and Prof. Gary Ayton • Help from Prof. Hans C. Andersen (Stanford) • I implemented the force-minimization approach * W. Noid, P. Liu, Y. Wang et al. J. Chem. Phys. 128, 244115 (2008).

22. Multiscale Coarse-Graining by Force Minimization Effective force: Central pairwise, linear approximation Each CG site: Residual: Multidimensional parabola Obtained from all-atom configurations

23. Force Minimization by Conjugate Gradient Method • Solving matrix directly Residual: Variational principle: • Or finding the minimal solution by conjugate gradient minimization with Ψ and gd • Subtract the Ewald Sum (long-range electrostatic) of point net charges • Match bonded and non-bonded interactions separately • Only one minimal solution! • Ψ can be used to determine the best CG scheme

24. Problems with MS-CG • EF-CG non-bonded effective forces Effective Force Coarse-Graining (EF-CG) Method* • Very limited transferability: temperature, surface, sequence of amino acids • Wrong pressure (density) without further constraint • Explicitly calculating pairwise atomic interactions between two groups • All-atom MD to get the ensemble of relative orientations * Y. Wang, W. Noid, P. Liu, G. A. Voth to be submitted.

25. Conclusions • CG methods enable faster simulations and longer effective simulation time • MS-CG method rebuilds structures accurately but has very limited transferability • MS-CG method can eliminate some atoms (e.g., implicit solvent) • EF-CG method has much better transferability by compromising a little accuracy of structures

26. MS-CG MD Study of Aggregation of Polyglutamines* • Polyglutamine aggregation is the clinic cause of 14 neural diseases, including Huntington’s, Alzheimer's, and Parkinson's diseases • All-atom simulations have a very slow dynamics that can not be adequately sampled • Water-free MS-CG model • CG MD simulations extend from nanoseconds to milliseconds CG MD results consistent with experiments: • Longer chain system exhibits stronger aggregation • Degrees of aggregation depend on concentration • Mechanism based on weak VDW interactions and fluctuation nature * Y. Wang, G. A. Voth to be submitted.

27. Some Applications of Ionic Liquids • Ionic liquid = Room temperature molten salt • Non-volatile • High viscosity Environment-friendly solvent for chemical reactions Lubricant Propellant

28. Multiscale Coarse-Graining of Ionic Liquids* • EMIM+/NO3- ionic liquid • 64 ion pairs, T = 400 K • Electrostatic and VDW interactions * Y. Wang, S. Izvekov, T. Yan, and G. Voth, J. Phys. Chem. B 110, 3564 (2006).

29. Site-site RDFs (T = 400K) • Good structures • No temperature transferability Satisfactory CG Structures of Ionic Liquids

30. Spatial Heterogeneity in Ionic Liquids* With longer cationic side chains: • Polar head groups and anions retain local structure due to electrostatic interactions • Nonpolar tail groups aggregate to form separate domains due to VDW interactions C1 C2 C4 C6 C8 * Y. Wang, G. A. Voth, J. Am. Chem. Soc. 127, 12192 (2005).

31. Define Heterogeneity order parameter (HOP) • For each site • Average over all sites to get <h> • Invariant with box size L Heterogeneity Order Parameter* • Larger HOP represents more heterogeneous configuration. • Quantifying degrees of heterogeneous distribution by a single value • Detecting aggregation • Monitoring self-assembly process * Y. Wang,G. A. VothJ. Phys. Chem. B110, 18601 (2006).

32. Thermal Stability of Tail Domain in Ionic Liquids* • Heat capacity plot shows a second order transition at T = 1200 K • Contradictory: HOP of instantaneous configurations do not show a transition at T = 1200 K? * Y. Wang, G. A. Voth, J. Phys. Chem. B 110, 18601 (2006).

33. Define Lattice HOP Mechanism • Instantaneous LHOPs at T = 1230 K • Divide simulation box into cells • In each cell the ensemble average of HOP is taken for all configurations • Heterogeneous tail domains have fixed positions at low T (solid-like structure) • Tail domains are more smeared with increasing T • Above Tc, instantaneous tail domains still form (liquid-like structure), but have a uniform ensemble average Tail Domain Diffusion in Ionic Liquids

34. Extendable EF-CG Models of Ionic Liquids* • CG force library • Extendibility, transferability, and manipulability • Extendable CG models correctly rebuild spatial heterogeneity features • CG RDFs do not change much for C12 from 512 (27,136) to 4096 ion pairs (217,088 atoms) • Proving spatial heterogeneity is truly nano-scale, not artificial effect of finite-size effect * Y. Wang, S. Feng, G. A. Voth J. Chem. Theor. Comp.5, 1091 (2009).

35. Disordering and Reordering of Ionic Liquids under an External Electric Field* From heterogeneous to homogeneous to nematic-like due to the effective screening of the external electric field to the internal electrostatic interactions. * Y. Wang J.Phys. Chem. B 113, 11058 (2009).

36. Conclusions • Spatial heterogeneity phenomenon was found in ionic liquids, attributed to the competition of electrostatic and VDW interactions • Solid-like tail domains in ionic liquids go through a second order melting-like transition and become liquid-like above Tc • EF-CG method was applied to build extendable and transferable CG models for ionic liquids, which is important for the systematic design of ionic liquids • Ionic liquid structure changes from spatial heterogeneous to homogeneous to nematic-like under an external electric field

37. Polymers for Gas-Separation Membranes UBE.com CO2 Capturer Air Dryer Air Mask • Environmental applications • Energy applications • Industrial applications • Military applications • …

38. Determining Crystalline Structure of Polymers Polybenzimidazole (PBI) • AMBER force field • Put one-unit molecules on lattice positions • Relax at P = 1 atm and T = 10 K • Measure lattice constants in relaxed configuration

39. Infinitely-Long Crystalline Polymers at T = 300 K X-Z Plane Y-Z Plane Polybenzimidazole (PBI) Poly[bis(isobutoxycarbonyl)benzimidazole] (PBI-Butyl) Kapton

40. CO2 and N2 inside PBI PBI + CO2 Very stiff PBI + N2 Sizes along Y are expanded. Gas molecules can hardly get in between the layers.

41. CO2 and N2 inside PBI-Butyl PBI-Butyl + CO2 Open up spaces PBI-Butyl + N2 No dimension sizes are changed. Gas molecules are free to diffuse between layers.

42. CO2 and N2 inside Kapton Kapton + CO2 Flexible Kapton + N2 Sizes along Z are expanded. Gas molecules change the crystal structure of Kapton.

43. Conclusions • PBI forms a very strong and closely packed crystalline structure. • CO2 and N2 can hardly diffuse in PBI crystal. • Crystal structure of PBI-Butyl is rigid, but the butyl side chains make the interlayer distances larger. • CO2 and N2 can freely diffuse between the layers. • Kapton crystal structure is also closely packed, but the interlayer coupling is weaker than in PBI. • CO2 and N2 can be accommodated between the layers which increases the interlayer distances. • CO2 and N2 behave similar in these three crystalline polymers.

44. Cracking of Crystalline PBI by Water (I) Water PBI Initial Final • Water molecules are attracted to PBI surface • Water molecules do not penetrate inside PBI • Water cluster suppresses the collective thermal vibration of PBI crystal

45. Cracking of Crystalline PBI by Water (II) 16 water molecules Initial Middle Final • Water molecules stick together by hydrogen bonds • PBI crystal structure change slightly

46. Cracking of Crystalline PBI by Water (III) Initial Final 160 water molecules • Water molecules form hydrogen bonding network • PBI crystal structure change drastically

47. Conclusions • To crack the crystal structure, PBI must have defects. • Strong binding of water molecules by hydrogen bonding network is possible to destroy local PBI crystal structures, thus to crack the crystal.

48. Fluctuation Theorems • Jarzynski’s equality: ensemble average over all nonequilibrium trajectories C. Jarzynski Phys. Rev. Lett. 78, 2690 (1997) • Crook’s theorem: involving nonequilibrium trajectories for both ways G. E. Crooks Phys. Rev. E 60, 2721 (1999) • Calculate free energy difference from fast nonequilibrium simulations. • Transiently absorb heat from environment.

49. 高级研究生课程 分子建模与模拟导论：2009年秋季 星期三下午15:20 – 17:00 S102教室