Surface Chemistry Affects the Efficacy of the Hydration Force between

3 hours ago - The hydration force is an important component of the energy landscape for particle-particle aggregation. Here we use molecular simulatio...
0 downloads 3 Views 1MB Size
Subscriber access provided by Kaohsiung Medical University

C: Surfaces, Interfaces, Porous Materials, and Catalysis

Surface Chemistry Affects the Efficacy of the Hydration Force between Two ZnO (101 #0) Surfaces Zhizhang Shen, Jaehun Chun, Kevin M. Rosso, and Christopher J. Mundy J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.8b02421 • Publication Date (Web): 18 May 2018 Downloaded from http://pubs.acs.org on May 18, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Surface Chemistry Affects the Efficacy of the Hydration Force between Two ZnO (1010) Surfaces

Zhizhang Shen,a* Jaehun Chun, a Kevin M. Rosso,a and Christopher J. Mundyab*

a

Physical Sciences Division, Physical & Computational Sciences Directorate, Pacific Northwest National Laboratory, Richland, Washington 99352, United States.

b

Department of Chemical Engineering, University of Washington, Seattle, Washington, United

States

*

Corresponding Author

Email: [email protected] [email protected]

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ABSTRACT: The hydration force is an important component of the energy landscape for particle-particle aggregation. Here we use molecular simulation to investigate the length-scale and the oscillatory nature of the hydration force that manifests between two hydrophilic surfaces that undergo hydroxylation in water. We identify how this force is modified due to the details of chemistry at the interface of ZnO (1010) in aqueous solution. Importantly, our research demonstrates that the size-dependent nature of oriented attachment does not originate from the hydration force. Rather, we find that the hydration force has qualitatively different outcomes based on the details of the surface chemistry occurring at structurally specific aqueous solidliquid interfaces.

ACS Paragon Plus Environment

Page 2 of 25

Page 3 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

INTRODUCTION Crystal growth via oriented attachment (OA) has been recognized in a wide range of materials, such as iron oxides1, zinc oxide2, titanium oxides3, iron (oxy)hydroxides4, gold5, platinum 6, and calcite7. Despite advances in nano-scale inter-particle force measurements8,9, theoretical understanding of the driving force is still not sufficiently comprehensive to predict OA outcomes. Conventional Derjaguin-Landau-Verwey-Overbeek (DLVO) theory commonly used for understanding colloidal particle-particle interactions fails to predict both the correct magnitude of the forces between nano-assemblies and the important angular dependence of OA10. Recent research suggests that the molecular details of solvent structure and ion hydration at the solid-liquid interface are needed to correct continuum-based frameworks such as DLVO, in order to provide a basis for co-alignment11. To this end, molecular dynamics (MD) simulations revealed the critical role of a non-uniform solution response for particle aggregation in aqueous environments5,9. According to these previous studies, adsorbed solvent creates higher energy barrier for coalescence when particles are misaligned at separations of 5 to 10 Å. MD simulations using reactive force fields showed the importance of the hydrogen bond network at separations < 5 Å, between surface hydroxyls for OA of metal oxides in humid environments12. This suggests that surface chemistry, i.e., the propensity to dissociate water, could dictate the occurrence of OA. However, to date, there is little research that provides a consistent picture between the work in Ref9 and Ref12. Hence the present study is directed at understanding the role of hydroxylation/surface chemistry on the interfacial solvent structure and the concomitant hydration force that nano-particles need to overcome to establish coalescence. In the present study the ZnO (1010)-water interface and the (1010) - (1010) interaction with intervening water was chosen as a model system. The ZnO (1010) surface is one of the most

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

commonly exposed and best understood surfaces of ZnO nano-particles. The surface chemistry of ZnO (1010) exposed up to several water layers is well characterized and both experimental and theoretical work find congruent results13,14. Although ZnO (0001) basal surfaces are a more likely candidate than (1010) for OA in experiments2,15, the study of latter avoids the complexity of reconstruction to eliminate the dipole on the former surface. Specifically, different reconstruction mechanisms have been revealed for Zn and O-terminated (0001) surfaces16–20, and a large discrepancy exists between experiments and calculations in terms of how much dissociation of adsorbed water occurs on the O-terminated surface21,22. Nevertheless, the study of the (1010) surface does not change the nature of the question we strive to resolve. Furthermore, no quantum mechanics-based method has been used to investigate the ZnO (1010)-bulk water interface, and such robust information is needed to establish both consistency and reliability of ReaxFF and other empirical force fields for accurately predicting key aspects of the aqueous interfacial chemistry, and thus the correct origins of the hydration force between two solid interfaces14,23. The term hydration force herein refers to water-mediated forces between two surfaces involving only the water-water and water-surface interactions24. In practice, the hydration force is obtained by subtracting the intrinsic particle-particle forces, such as the dispersion force, from the total force between particles in the aqueous environment. The importance of the hydration force for hydrophilic surfaces is that the repulsive hydration force dominates at short separations24. For hydrophobic surfaces, the attractive hydration force prevails when the two approach each other, in accordance with the hydrophobic effect25. For hydrophilic hard surfaces (e.g., mica), the hydration force is usually composed of an oscillatory force related to the ordering of water molecules into layers, due to the presence of a surface creating an enveloping

ACS Paragon Plus Environment

Page 4 of 25

Page 5 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

repulsive force with a monotonic decay26. However, a recent molecular simulation revealed both attractive and repulsive hydration force between hydrophilic surfaces, highlighting the interplay between direct surface interactions and the water binding affinity of the surfaces27. In this article, we first report the surface hydroxylation and water structure at the ZnO (1010)-bulk water interface, calculated from quantum density functional theory molecular dynamics (DFT-MD). To see the effect of hydroxylation on the surface water structure and the assembly process, we also include the results from a modified DFT-MD (no hydroxylation) and from classical MD using both reactive force fields (ReaxFF)28 and non-reactive force fields9. Because of the computational resources needed for the system sizes required to compute the potential of mean force (PMF) between two surfaces using DFT-MD, we resort to ReaxFF and non-reactive force fields to calculate the PMFs in a slab-particle setup (Fig. 1(a)). Understanding the correlation between the response of water to the ZnO (1010) surface, the hydration force, and the changes due to the chemistry, namely hydroxylation, are important outcomes of this study. Since OA occurs more often via end-to-end attachment with smaller surface areas5, we also attempt to address the size/shape effect on the hydration force by simulating three different systems in this study, namely, the slab-particle, particle-particle, and slab-slab systems (Fig. 1).

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 1: Three simulations setups used for PMF calculations in this work: (a) slab-particle system, (b) particle-particle system, and (c) slab-slab system. THEORETICAL METHODS Non-reactive Force Fields MD Simulation Setup. MD simulations at room temperature were performed by using DL_POLY Classic29 in NVT (constant number of particles, constant volume, and constant temperature) ensemble. The Nosé-Hoover thermostat30 was used with a relaxation time of 0.5 ps. The equations of motion were integrated using the Verlet-leapfrog algorithm with a time step of 2 fs. Real-space cutoffs of 9.35 Å, 10 Å, and 10 Å were used in the Ewald sum for the slab-slab, slab-particle, and particle-particle systems, respectively. The same cutoffs were used for the short-range interactions in respective system. The slab-slab system comprises two of 6 double-layer slabs with x-y dimensions of 19.625 × 20.552 Å2. In the slabparticle system, the 6 double-layer slab substrate, with x-y dimensions of 52.333 × 51.381 Å2, interacted with (100) surface of ZnO particle on the top. This ZnO particle has a shape of hexagonal prism with six (100) surfaces and two (001) surfaces. The “contact area” between the substrate and the particle is about 12.66 × 16.02 Å2. This same particle was used in the particleparticle system. The SPC/E model was used to simulate water molecules31. The force fields for ZnO and interactions between ZnO and SPC/E water were obtained from reference9. GCMC Simulations. The amount of water molecules in between two ZnO slabs at each separation at 300 K was determined by using configurational bias grand canonical Monte Carlo (constant µVT) simulations32, where µ is the water chemical potential. The simulations were performed by using the Towhee code33. The system consists of a 10 double-layer of ZnO (with the same x-y dimensions as slab-slab MD simulations), whose top and bottom (100) surfaces were separated by distances ranged from 5 Å to 40 Å. The coordinates of Zn and O atoms were kept fixed since trivial surface reconstructions were observed in MD simulations, while the water

ACS Paragon Plus Environment

Page 6 of 25

Page 7 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

molecules could translate and rotate, and equilibrate with a SPC/E water reservoir with a chemical potential of - 49.0 kJ/mol at room temperature34. The following probabilities were used for GCMC moves: 40% configurational-bias insertion or deletion, 15% intrabox configurationalbias molecule transfer, 23% center-of-mass translation, and 22% rotation about the center-ofmass. For separations that are less than or equal to 20 Å, simulations were carried out for a total of 3 × 107 moves, of which the last 107 moves were used for obtaining the average number of water molecules. In order to promote the sampling efficiency for separations larger than 20 Å, each system was first equilibrated with a chemical potential that is 10% higher than that of SPC/E water for 4 × 106 steps and then with the correct chemical potential for 4 × 107 steps. ReaxFF MD Simulations Setup. The same slab-particle system was used in ReaxFF MD calculations except that the x-y dimensions were changed to 52.74 × 53.10 Å2 based on the cell parameters optimized by ReaxFF35. LAMMPS code36 with the extended package USERREAXC37 was used for all calculations. A time step of 0.25 fs was used. The Nosé-Hoover thermostat was used to kept temperature constant in NVT ensemble with a relaxation time of 0.5 ps. Potential of Mean Force (PMF) Calculations. In slab-slab and slab-particle systems, the center-of-masses (COMs) of both objects along z-direction were constrained but both objects were allowed to translate on x-y plane. In particle-particle system, the COMs along all three dimensions were constrained. In all cases, the rotation of an object about the COM was removed. The average constraining normal fz of a whole slab or particle was used to calculate the potential of mean force as a function of separation between two objects, slab or particle, by the following



equation: Ar = 〈f 〉dr ≈ 〈f 〉dr, where r is the distance between the geometric center of 

the two facing surface layers of the two objects along the z-direction, and z0 was taken as 20 Å

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

for slab-slab system or 15 Å for the other two systems in the calculations using non-reactive force fields. In each constrained MD calculation by using non-reactive force fields, the system was run for 3ns and the constraining force was collected for the last 1ns. The hydration force fz,hdy in each case was obtained by subtracting the object-object force fz,o-o in vacuum from the total force fz: fz,hdy = fz - fz,o-o. The object-object force fz,o-o was acquired by MD simulations of each system with water molecules being removed, i.e., in vacuum. In the PMF calculations using ReaxFF, the initial configuration at each separation was taken from the corresponding structure pre-equilibrated for 1ns using non-reactive force fields. z0 was taken as 10 Å in the ReaxFF PMF calculations. Each system was run for 400 ps and the last 200 ps was used for collecting forces. DFT MD Simulations. Born-Oppenheimer MD simulations were performed by using the program CP2K/Quickstep38. The initial structure for DFT-MD was obtained from a 1 ns long classical MD simulation using non-reactive force fields where the surface slab is composed of six Zn-O double layers with a surface area of 19.625 × 20.553 Å2. The two surfaces were separated by a water box of 25 Å thickness. The simulations were run in canonical ensemble (NVT) using a Nosé-Hoover thermostat with a target temperature of 300 K. We used revised functional of Perdew-Burke-Ernzerhof (revPBE39) with D3 (Grimme’s third generation of dispersion40). This exchange-correlation functional was utilized because it has been shown to correctly describe the properties of water at ambient conditions41 while the water adsorption energy and surface chemistry of ZnO (1010) are not sensitive to which type of dispersion is used and whether a fraction of exact exchange is included14. We used a cutoff of 400 Ry for the plan wave expansion and double- valence polarized basis sets optimized to the condensed phase. Twelve valence electrons were used for Zn (3d104s2), six for O (2s22p4), and one for H. Only the

ACS Paragon Plus Environment

Page 8 of 25

Page 9 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

gamma point in the Brillouin zone was used. We carried out a total of 42 ps (20 ps of equilibration run and 22 ps of production run) and 30 ps (15 ps of equilibration run and 15 ps of production run) for unconstrained and constrained MD runs, respectively, with a time step of 0.5 fs. In the constrained MD run, we controlled the hydroxylation percentage at zero by applying a quadratic restraining potential to a collective variable (CV). This CV can be expressed as a 

Gaussian function of coordination number as follows: POP = ∑#$%& ' √ e



 !   "

. Here n) =

0, σ = 0.5 and cn is the coordination number variable that is a simple function of the distance between surface oxygen and hydrogen atoms defined by cn = ∑/1 − r⁄r) .. ⁄1 − r⁄r) .0 , where r) = 1.98, nn = 12, nd = 24. By using the parameters above, the population CV representing zero hydroxylation is 0.8 for each surface O site and 19.2 for each surface (24 sites in total). The applied quadratic restraining potential is in the form of U = kPOP − cv , where k = 10 Hartree, and cv = 19.2. RESULTS 50) Surface. To investigate the change in the solvent response Water Structure above (104 due to hydroxylation and then to contrast the result from force fields to the DFT-based approaches, we performed four sets of simulations of the ZnO (1010)-water interface (Fig. 2(a)) using ReaxFF28, non-reactive force fields9, DFT-MD, and modified DFT-MD in which the degree of surface hydroxylation is restrained at 0%. The resulting density profiles are shown in Fig. 2(b) revealing differences between the choices of interaction potential. At least three pronounced water layers exist in all models. However, it is clear from Fig. 2(b) that the structure of first water layer is qualitatively changed by surface hydroxylation. First, the bimodal distribution of the first water layer in both reactive surface models (ReaxFF and DFT-MD) is

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

caused by surface water dissociation. The unreactive surface gives rise to a hydrophobic region (or density lower than that of the bulk water) at 2 Å that is present for both unreactive simulations. It is important to note that the solution response resulting from the unreactive FF depicts similar peak positions and an overall profile shape in good qualitative agreement with the modified DFT-MD counterpart. In both reactive models, three similar types of proton hopping occur on the surface (Fig. 2(c), (d), (e)). The free energy surface for proton transfer to and from the surface and transfer between adsorbed water and OHs can be estimated from DFT-MD trajectories. The free energy is related to the probability distribution of two distances (do-o and doh-oh) by ∆7 = −89 : log > where do-o is the distance between two oxygens as H donor and acceptor, and doh-oh is the difference in the distances between H and the two oxygens. From the proton transfer statistics in DFT-MD, free energy barriers for hopping to and from the surface and for hopping between adsorbed OHs and H2Os are 0.04 and 0.06 eV, respectively. These values are lower than the corresponding values in a water film adsorption system (0.07 and 0.1 eV, respectively14). The hydroxylation percentage, or the percentage of surface O atoms bound to H atoms, in the DFT-MD is 62.5%, whereas the ReaxFF simulations produce 80~85% surface hydroxylation. The ReaxFF curve has an additional water layer that cannot be distinguished from the shoulder peak of the first layer, in contrast to the DFT representation of interaction (gray region and the inset of Fig. 2(b)). Overall, all four models predict qualitatively similar water response above the surface, with the nonreactive FF and ReaxFF producing the most ordered and disordered density profiles, respectively.

ACS Paragon Plus Environment

Page 10 of 25

Page 11 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 2. (a) A schematic of the water structuring at the surface of ZnO (1010). Zinc is shown in gray, oxygen in red, and hydrogen in white. (b) Water density profiles on ZnO (1010) calculated from DFT MD (black), modified DFT MD (green), ReaxFF (red), and the non-reactive FF (blue). The curve for the non-reactive FF is shifted left by 0.3 Å to better resolve the difference between representations of the interaction. The region of the second water layer is shaded in gray. The decomposition of O peaks of the first two water layers for ReaxFF is shown in the inset. Note that the oxygen distributions of the 1st and 2nd water layers overlap. The second water layer O density is taken as the difference between the overall O density and the first water layer O density. (c) For the DFT-MD and ReaxFF, protons hop back-andforth between molecularly adsorbed water and surface O atoms. (d) Proton transfer between the neighboring adsorbed H2Os and OHs. (e) Grotthus type proton transfer.

50) Surfaces. How solvent response to two bound surfaces PMFs between Two (104 influences the strength and length scale of the hydration force was then investigated, by using non-reactive force fields and the ReaxFF, which represent the two end-members of the surface

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

chemistry spectrum (0% and 85% hydroxylation). Focusing on the slab-particle model first, PMF curves for both the total force and hydration force using the non-reactive force fields are shown in Figs 3(a) and (b). As seen from Fig. 3(a), the hydration force dominates nearly the entire range of separation explored over the intrinsic slab-particle forces. It is important to note that the hydration force predominance is due to the inherent short-range nature of the empirical nonbonded interaction that does not contain a high fidelity long-range interaction (~ O(1) nm), which would be described by a modified Hamaker function as discussed in Ref11. Our assumption is that the aqueous response to solid interfaces as described by the water—surface interaction potentials used herein does provide a good description of hydration and is consistent with assumptions of other related studies5,9. Three attractive hydration force maxima (at 5.75 Å, 8 Å, and 10 Å, respectively) are also revealed from Fig. 3(a). These maxima exist between two surfaces through the hydrophobic water structure shown in Figs. 3(d) and (f). The water structure between two unbound surfaces is also included as a reference (Fig. 3(c)). Similarly, repulsive hydration forces occur through elevated water density regions, or hydrophilic regions occurring when the particle density profiles overlap (Fig. 3(e)). In the total PMF (Fig. 3(b)), a deep energy well of -0.65 eV occurs at ~ 7.6 Å when two water layers reside on each surface of the two particles. If the particles are moved closer a seemingly insurmountable free energy barrier of ~ 2.5 eV is encountered, making contact a rare event (gray region in Fig. 3(b)). The energy barrier for moving from the well at ~10 Å (three water layers reside on each surface) to the deep well at ~ 7.6 Å is about 0.16 eV (cyan region in Fig. 3(b)).

ACS Paragon Plus Environment

Page 12 of 25

Page 13 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 3. (a,b) Total (in red) and hydration (in black) force curves for slab-particle system, calculated from the non-reactive FF and the corresponding PMFs. (c, d, e, f) Water O density profiles at separations of 40, 10, 9, and 8 Å, respectively. The x-axis in each density profile is the separation between two

surfaces in Å with an arbitrary zero distance and the y-axis is the O density in g/cm3.

Figure 4. (a,b) Total (in black) and hydration (in red) force curves for the same system using ReaxFF and the corresponding PMFs, respectively. The inset in (b) shows the enlarged PMF curve at separations from 7.5 to 10 Å.

Similar PMFs are accessible using ReaxFF for the identical slab-particle system, as shown in Fig. 4(a) and (b). Both the sign and magnitude of the hydration force correlate well with the

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

water structure on the surface. Based on the total PMF curve (Fig. 4(b)), it is almost barrier-less for the particles to overcome the second water layer on each surface of the two slabs to move to the “locked-up” position, where only one water layer remains on each surface (gray region in Fig. 4(b)). Two energy barriers with similar values (0.10~0.20 eV) exist for moving from the separation where three water layers reside on each surface to the one where only two water layers are retained (cyan region in Fig. 4(b)). Size Dependence of Hydration Force. Comparison of these results from the slab-particle model to those using a semi-infinite slab-slab approach versus a particle-particle approach was facilitated using the efficient non-reactive force fields. We can then test these various simulation systems using intuition for how the hydration forces differ for both particle size and shape. For the PMF between two semi-infinite slabs, the number density of water for each slab-slab separation was determined from configurational bias GCMC calculations. In order to compare the effect of the interparticle force on system geometry the calculated thermodynamic forces for the three systems are normalized to the surface area (i.e., pressure) and plotted in Fig. 5. In the slab-slab and particle-particle systems, the hydration force remains the dominant force between two surfaces. The hydration pressure curves are consistent and nearly identical at separations > 7 Å and suggests that the hydration force scales linearly with the surface area and that all simulation geometries yield consistent interactions.

ACS Paragon Plus Environment

Page 14 of 25

Page 15 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 5. Hydration pressure as a function of separation for slab-slab (black), slab-particle (red), and particle-particle (blue) systems. Differences occur in the magnitude of the peak at ~6 Å that might be caused by a GCMC sampling issue at close separations. The schematics of these three systems are shown at the right. The dotted lines schematically delineate the periodic boundary conditions used in each system.

DISCUSSION In principle, self-assembly/aggregation of nanoparticles (e.g., OA) is a result of a delicate balance between colloidal forces between the ZnO nanoparticles (i.e., hydration, electrostatic, and dispersion forces) and hydrodynamic forces42. Because standard empirical interaction potentials are short-ranged, dispersion forces between two surfaces in solution are typically modeled using continuum-based Lifshitz formulation that rely on the macroscopic property, namely the frequency-dependent dielectric function of ZnO. Based on the aforementioned formulation, the dispersion interaction between ZnO surfaces would scale as ~ 0.1 eV from a

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

non-retarded Hamaker function that is in a similar range to the well-studied mica system43. However, a recent study44 determined that for nanoparticles in solution, the effect of solution response to the solid interfaces can increase the dispersion interaction and thus suggests that the Hamaker function can be increased by one order of magnitude higher for separations of interest in our study. Therefore, the dispersion interaction, ~ O(1) eV, becomes comparable to the pure hydration interaction predicted in this study, underscoring the importance of needing different theoretical frameworks to elucidate both interactions to ensure a deep understanding of the forces contributing to OA. In addition to dispersion and hydration interactions, an electrostatic interaction can play a role in the assembly process. Using continuum electrostatics, the effective surface potential estimated from the measured zeta potential (i.e., ~ 40 mV at pH=6-745) is not expected to be large enough (i.e., mica is known to be ~ -150 mV of the effective surface potential46) to play a significant role at the length scale of interest (i.e., O(1) nm). Therefore, it is the dispersion interaction between ZnO nanoparticles that is a dominating attractive macroscopic interaction. Indeed, our study reveals that the assembly of ZnO nanoparticles will be dictated by a balance between dispersion and hydration interactions. Furthermore, as shown in Figure 5, no clear indication of a strong size dependence that alters the qualitative nature of the computed hydration interaction is evident, suggesting that the size-dependent nature of OA (i.e., that OA is the domain of small particles only and not macroscopic crystals) originates instead from other physiochemical nature. Such nature would include the role of surface roughness at close separations, diffusive transport associated with expulsion of intervening water, or the relative importance of solvent structuring effects on dispersion and/or electrostatic forces as to crystal sizes.

ACS Paragon Plus Environment

Page 16 of 25

Page 17 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Experimental OA studies of metal oxides in humid environments can provide useful information regarding the final stage of aggregation8,12. Based on different theoretical frameworks, namely MD or DLVO, both van der Waals torque and hydrogen bonding can regulate the alignment at this close proximity (< 5 Å). The present study also connects to the complementary process of the necessary removal of water that is needed to achieve OA. The “locked-up” configuration at 5~6 Å between ZnO nanoparticles where the adsorbed water/hydroxyls of one nanoparticle hydrogen-bond to the surface O or hydroxyls of another aligned nanoparticle roughly resemble the two aligned nanoparticles connected by hydrogen bonds in a previous ReaxFF MD work applied to TiO212. At larger separations in aqueous environments, both the reactive and non-reactive representations of interaction reveal the existence of a metastable solvent-separated regime at ~10 Å (three water layers on each surface) and similar energy barrier height (0.1~0.2 eV) for moving from this solvent-separated regime to the solvent-separated regime at 7~8 Å (two water layers on each surface). However, there is a strong dependence on the representation of interactions regarding the predictions moving from ~7 to 8 Å to the “locked-up” position. Specifically, the non-reactive representations of interaction suggest that particles will mostly reside at the solvent separated regime at ~8 Å due to an energy barrier of 2.5 eV, whereas the ReaxFF result in no free-energy barrier to reach the “locked-up” position. This is dictated by surface hydroxylation resulting in a less ordered solvent response and leading to a barrier-less approach to the penultimate stage for nanoassembly. We believe that the observed chemistry of the ZnO (1010) surface is likely best described as a case between the limiting cases of both the non-reactive and ReaxFF representations of interaction. This difference in hydration free energy will likely influence the level of perfection of OA process where particles can align themselves in a stable solvent-separated regime that allow the

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

sufficient sampling of rotational motions. The hydrophilic case studied with a non-reactive FF is likely to display a solvent-separated regime between particles (at ~8 Å separation) and allow particles to sample rotational configurations. In contrast, a reactive interaction potential that describes chemistry at the surface such as ReaxFF will not produce a solvent-separated regime due to attractive nature of both hydration and dispersion interactions and will facilitate hydrophobic collapse as the dominant pathway to association. Our study demonstrates that the surface chemistry is critical to tuning the detailed nature of the resulting hydration force that in turn controls the pathway of the assembly process. To the extent that the current models used herein capture the correct surface phenomena (chemistry) that ultimately tunes the surface hydrophobicity requires additional research. Indeed, the tuning of surface hydrophobicity and hydrophilicity has been shown in model systems to produce dramatically different outcomes assembly that corroborate our current findings in the more complex system of ZnO47–49. This work provides a step to isolate the role of the hydration force in the nanoparticle assembly and provide both the necessary simulation protocol and the details of the sensitivity of the hydration force due to the representation of interaction. Our work suggests that DFT-based PMF calculations between ZnO nanoparticles in solution, which are computationally intensive, are an imperative next step in order to accurately describe the strength of hydration force at separations < 10 Å. ACKNOWLEDGEMENT MD and DFT simulations and manuscript preparation were supported by the MS3 (Materials Synthesis and Simulation Across Scales) Initiative at Pacific Northwest National Laboratory (PNNL), a multi-program national laboratory operated by Battelle under Contract DE-AC0576RL01830 for the U.S. Department of Energy. Supporting MD/DFT simulations, data analysis,

ACS Paragon Plus Environment

Page 18 of 25

Page 19 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

and manuscript preparation were supported by the U.S. Department of Energy (DOE) office of Basic Energy Sciences (BES), Division of Chemical Sciences, Geosciences, and Biosciences. Theoretical interpretation for macroscopic forces/interactions was supported by the U.S. DOE BES, Division of Materials Sciences and Engineering. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. The research was partially performed using PNNL Institutional Computing. The authors would like to thank Evgenii Fetisov at the University of Minnesota for fruitful discussion regarding the GCMC simulations. REFERENCES (1)

Almeida, T. P.; Fay, M. W.; Zhu, Y.; Brown, P. D. In Situ TEM Investigation of βFeOOH and α-Fe2O3 Nanorods. Phys. E Low-dimensional Syst. Nanostructures 2012, 44, 1058–1061.

(2)

Pacholski, C.; Kornowski, A.; Weller, H. Self-Assembly of ZnO: From Nanodots to Nanorods. Angew. Chemie - Int. Ed. 2002, 41, 1188–1191.

(3)

Penn, R. L. Imperfect Oriented Attachment: Dislocation Generation in Defect-Free Nanocrystals. Science 1998, 281, 969–971.

(4)

Li, D.; Nielsen, M. H. M. H.; Lee, J. R. I. J. R. I.; Frandsen, C.; Banfield, J. F.; De Yoreo, J. J. Direction-Specific Interactions. Science 2012, 336, 1014–1018.

(5)

Welch, D. A.; Woehl, T. J.; Park, C.; Faller, R.; Evans, J. E.; Browning, N. D. Understanding the Role of Solvation Forces on the Preferential Attachment of Nanoparticles in Liquid. ACS Nano 2016, 10, 181–187.

(6)

Zheng, H.; Smith, R. K.; Jun, Y. W.; Kisielowski, C.; Dahmen, U.; Alivisatos, A. P.

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Observation of Single Colloidal Platinum Nanocrystal Growth Trajectories. Science 2009, 324, 1309–1312. (7)

Takasaki, M.; Oaki, Y.; Imai, H. Oriented Attachment of Calcite Nanocrystals: Formation of Single-Crystalline Configurations as 3D Bundles via Lateral Stacking of 1D Chains. Langmuir 2017, 33, 1516–1520.

(8)

Zhang, X.; He, Y.; Sushko, M. L.; Liu, J.; Luo, L.; De Yoreo, J. J.; Mao, S. X.; Wang, C.; Rosso, K. M. Direction-Specific van Der Waals Attraction between Rutile TiO 2 Nanocrystals. Science 2017, 356, 434–437.

(9)

Zhang, X.; Shen, Z.; Liu, J.; Kerisit, S. N.; Bowden, M. E.; Sushko, M. L.; De Yoreo, J. J.; Rosso, K. M. Direction-Specific Interaction Forces Underlying Zinc Oxide Crystal Growth by Oriented Attachment. Nat. Commun. 2017, 8, 835.

(10)

De Yoreo, J. J.; Gilbert, P. U. P. A.; Sommerdijk, N. A. J. M.; Penn, R. L.; Whitelam, S.; Joester, D.; Zhang, H.; Rimer, J. D.; Navrotsky, A.; Banfield, J. F.; et al. Crystallization by Particle Attachment in Synthetic, Biogenic, and Geologic Environments. Science 2015, 349, aaa6760.

(11)

Li, D.; Chun, J.; Xiao, D.; Zhou, W.; Cai, H.; Zhang, L.; Rosso, K. M.; Mundy, C. J.; Schenter, G. K.; De Yoreo, J. J. Trends in Mica–Mica Adhesion Reflect the Influence of Molecular Details on Long-Range Dispersion Forces Underlying Aggregation and Coalignment. Proc. Natl. Acad. Sci. 2017, 114, 7537–7542.

(12)

Raju, M.; Van Duin, A. C. T.; Fichthorn, K. A. Mechanisms of Oriented Attachment of TiO2 Nanocrystals in Vacuum and Humid Environments: Reactive Molecular Dynamics. Nano Lett. 2014, 14, 1836–1842.

(13)

Meyer, B.; Rabaa, H.; Marx, D. Water Adsorption on ZnO(1010): From Single Molecules

ACS Paragon Plus Environment

Page 20 of 25

Page 21 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

to Partially Dissociated Monolayers. Phys. Chem. Chem. Phys. 2006, 8, 1513. (14)

Tocci, G.; Michaelides, A. Solvent-Induced Proton Hopping at a Water − Oxide Interface. J. Phys. Chem. Lett. 2014, 5, 474–480.

(15)

Fan, B.; Zhang, Y.; Yan, R.; Fan, J. Multistage Growth of Monocrystalline ZnO Nanowires and Twin-Nanorods: Oriented Attachment and Role of the Spontaneous Polarization Force. CrystEngComm 2016, 18, 6492–6501.

(16)

Dulub, O.; Diebold, U.; Kresse, G. Novel Stabilization Mechanism on Polar Surfaces: ZnO(0001)-Zn. Phys. Rev. Lett. 2003, 90, 16102.

(17)

Kresse, G.; Dulub, O.; Diebold, U. Competing Stabilization Mechanism for the Polar ZnO(0001)-Zn Surface. Phys. Rev. B 2003, 68, 1–15.

(18)

Önsten, A.; Stoltz, D.; Palmgren, P.; Yu, S.; Göthelid, M.; Karlsson, U. O. Water Adsorption on ZnO(0001): Transition from Triangular Surface Structures to a Disordered Hydroxyl Terminated Phase. J. Phys. Chem. C 2010, 114, 11157–11161.

(19)

Lauritsen, J. V.; Porsgaard, S.; Rasmussen, M. K.; Jensen, M. C. R.; Bechstein, R.; Meinander, K.; Clausen, B. S.; Helveg, S.; Wahl, R.; Kresse, G.; et al. Stabilization Principles for Polar Surfaces of ZnO. ACS Nano 2011, 5, 5987–5994.

(20)

Kunat, M.; Girol, S. G.; Burghaus, U.; Wöll, C. The Interaction of Water with the Oxygen-Terminated, Polar Surface of ZnO. J. Phys. Chem. B 2003, 107, 14350–14356.

(21)

Wahl, R.; Lauritsen, J. V.; Besenbacher, F.; Kresse, G. Stabilization Mechanism for the Polar ZnO(0001)-O Surface. Phys. Rev. B - Condens. Matter Mater. Phys. 2013, 87, 1–12.

(22)

Schiek, M.; Al-Shamery, K.; Kunat, M.; Traeger, F.; Wöll, C. Water on the Hydroxylated H-(1×1) O-ZnO (0001) Surface. Phys. Chem. Chem. Phys. 2006, 8, 1505–1512.

(23)

Große Holthaus, S.; Köppen, S.; Frauenheim, T.; Colombi Ciacchi, L.; Holthaus, S.; Ko,

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

S.; Frauenheim, T.; Ciacchi, L. C.; Große Holthaus, S.; Köppen, S.; et al. Atomistic Simulations of the ZnO(1210)/Water Interface: A Comparison between First-Principles, Tight-Binding, and Empirical Methods. J. Chem. Theory Comput. 2012, 8, 4517–4526. (24)

Parsegian, V. A.; Zemb, T. Hydration Forces: Observations, Explanations, Expectations, Questions. Curr. Opin. Colloid Interface Sci. 2011, 16, 618–624.

(25)

Meyer, E. E.; Rosenberg, K. J.; Israelachvili, J. Recent Progress in Understanding Hydrophobic Interactions. Proc. Natl. Acad. Sci. U. S. A. 2006, 103, 15739–15746.

(26)

Kilpatrick, J. I.; Loh, S. H.; Jarvis, S. P. Directly Probing the Effects of Ions on Hydration Forces at Interfaces. J. Am. Chem. Soc. 2013, 135, 2628–2634.

(27)

Kanduč, M.; Schneck, E.; Netz, R. R. Attraction between Hydrated Hydrophilic Surfaces. Chem. Phys. Lett. 2014, 610, 375–380.

(28)

Raymand, D.; van Duin, a. C TRaymand, D., van Duin, a. C. T., Baudin, M., & Hermansson, K. A Reactive Force Field (ReaxFF) for Zinc Oxide. Surf. Sci. 2008, 602, 1020–1031.

(29)

Smith, W.; Forester, T. R. DL_POLY_2.0: A General-Purpose Parallel Molecular Dynamics Simulation Package. J. Mol. Graph. 1996, 14, 136–141.

(30)

Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A 1985, 31, 1695–1697.

(31)

Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 6269–6271.

(32)

Mooij, G. C. A. M.; Frenkel, D.; Smit, B. Direct Simulation of Phase Equilibria of Chain Molecules. J. Phys. Condens. Matter 1992, 4, L255-L259.

(33)

Martin, M. G. MCCCS Towhee: A Tool for Monte Carlo Molecular Simulation. Mol.

ACS Paragon Plus Environment

Page 22 of 25

Page 23 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Simul. 2013, 39, 1212–1222. (34)

Boţan, A.; Rotenberg, B.; Marry, V.; Turq, P.; Noetinger, B. Hydrodynamics in Clay Nanopores. J. Phys. Chem. C 2011, 115, 16109–16115.

(35)

Raymand, D.; van Duin, A. C. T.; Spångberg, D.; Goddard, W. A.; Hermansson, K. Water Adsorption on Stepped ZnO Surfaces from MD Simulation. Surf. Sci. 2010, 604, 741–752.

(36)

Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1–42.

(37)

Aktulga, H. M.; Campus, M.; Fogarty, J. C.; Grama, A. Y. Parallel Reactive Molecular Dynamics : Numerical Methods and Algorithmic Techniques. Parallel Comput. 2012, 38, 245–259.

(38)

Vandevondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. Quickstep: Fast and Accurate Density Functional Calculations Using a Mixed Gaussian and Plane Waves Approach. Comput. Phys. Commun. 2005, 167, 103–128.

(39)

Zhang, Y.; Yang, W. Comment on “Generalized Gradient Approximation Made Simple.” Phys. Rev. Lett. 1998, 80, 890–890.

(40)

Grimme, S.; Stephan, E.; Goerigk, L. Effect of the Damping Function in Dispersion Corrected Density Functional Theory. J. Comput. Chem. 2011, 32, 1456–1465.

(41)

Galib, M.; Duignan, T. T.; Misteli, Y.; Baer, M. D.; Schenter, G. K.; Hutter, J.; Mundy, C. J. Mass Density Fluctuations in Quantum and Classical Descriptions of Liquid Water. J. Chem. Phys. 2017, 146, 244501.

(42)

Pednekar, S.; Chun, J.; Morris, J. F. Simulation of Shear Thickening in Attractive Colloidal Suspensions. Soft Matter 2017, 13, 1773–1779.

(43) Bergström, L. Hamaker Constants of Inorganic Materials. Adv. Colloid Interface Sci. 1997,

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

70, 125–169. (44)

Chun, J.; Mundy, C. J.; Schenter, G. K. The Role of Solvent Heterogeneity in Determining the Dispersion Interaction between Nanoassemblies. J. Phys. Chem. B 2015, 119, 5873– 5881.

(45)

Luo, M.; Shen, C.; Feltis, B. N.; Martin, L. L.; Hughes, A. E.; Wright, P. F. A.; Turney, T. W. Reducing ZnO Nanoparticle Cytotoxicity by Surface Modification. Nanoscale 2014, 6, 5791–5798.

(46)

Prakash, A.; Pfaendtner, J.; Chun, J.; Mundy, C. J. Quantifying the Molecular-Scale Aqueous Response to the Mica Surface. J. Phys. Chem. C 2017, 121, 18496–18504.

(47)

Eun, C.; Berkowitz, M. L. Fluctuations in Number of Water Molecules Confined between Nanoparticles. J. Phys. Chem. B 2010, 114, 13410–13414.

(48)

Acharya, H.; Vembanur, S.; Jamadagni, S. N.; Garde, S. Mapping Hydrophobicity at the Nanoscale. Faraday Discuss. 2010, 146, 353–365.

(49)

Patel, A. J.; Varilly, P.; Chandler, D. Fluctuation of Water near Extended Hydrophobic and Hydrophilic Surfaces. J. Phys. Chem. B 2010, 114, 1632–1637.

ACS Paragon Plus Environment

Page 24 of 25

Page 25 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

TOC Graphic

ACS Paragon Plus Environment