Article pubs.acs.org/JPCC
Reactive Force Field Modeling of Zinc Oxide Nanoparticle Formation Craig J. Tainter and George C. Schatz* Department of Chemistry, Northwestern University, Evanston, Illinois 60208 United States S Supporting Information *
ABSTRACT: A ReaxFF reactive force field is developed and used for molecular dynamics studies of reactions that occur when diethyl zinc is exposed to a graphene surface that has been functionalized with epoxide groups. From past experiments, it is known that these conditions lead to zinc oxide nanoparticle formation. Molecular dynamics simulations are used to provide atom-level detail into the nanoparticle formation process, including the mechanisms whereby oxygen is abstracted from the graphene surface, thus enabling condensation reactions in which multiple zinc-containing species form zinc oxide fragments and ultimately nanoparticles. Structural properties of the nanoparticles show nonstoichiometric zinc oxide structures with average coordination numbers of 3.6 around each zinc. Time-dependent density functional theory calculations show that the absorption spectra of these clusters are red-shifted by a few tenths of an electronvolt compared to that of a wurtzite crystal structure, representing transitions from oxygen 2p’s near the cluster surface to zinc 4s’s in the interior.
■
INTRODUCTION Metal oxide nanomaterials have been the focal point of many recent scientific studies because of their mechanical, electrical, and optical properties. Zinc oxide is one such semiconducting material that has received much attention because of its large band gap (3.36 eV) and large exciton binding energy (60 meV) as well as its good photoelectric and piezoelectric properties. These properties have led to many applications for zinc oxide (ZnO) including sensors,1−4 catalysts,5,6 and optoelectronic devices.7−9 The bulk crystal structure for zinc oxide is known to adopt the wurtzite structure; however, experiments have shown that high pressures can induce a transition to a cubic rocksalt structure.10−12 Furthermore, different nanostructures can be constructed by employing different experimental conditions and/or methods. Zinc oxide nanoparticles in the 2−7 nm size range can be formed using solution-based methods such as the addition of lithium hydroxide to an ethanolic zinc acetate solution.13 The particle size can be controlled by varying the temperature and concentration of different species during nanoparticle growth. Laser ablation of zinc metal in different surfactant solutions has also been used to grow ZnO nanoparticles.14 Here, roughly round nanoparticles in the 10−30 nm size range were formed in cationic, amphoteric, and nonionic surfactant solutions as well as in deionized water but were not seen in anionic surfactant solutions. Vapor transport techniques are also widely used methods to synthesize ZnO nanoparticles. A variety of nanostructures can be selectively formed from the thermal evaporation of ZnO powders at controlled growth conditions.15 In a 2005 review, Fan and Lu report many other vapor transport methods for producing ZnO nanoparticles as well as patterned growth using photolithography.16 © XXXX American Chemical Society
In this study we will focus on the recent article by Johns et al. which uses an atomic layer deposition (ALD) technique to form ZnO nanoparticles.17 In the Johns work, molecular oxygen was thermally cracked using a tungsten filament and the resulting atomic oxygen was deposited onto a graphene sheet. This activates graphene by functionalizing the surface with epoxide groups. Diethyl zinc (henceforth denoted EZE) was subsequently deposited onto this surface after which atomic force microscopy (AFM) images show nanoparticle (presumably ZnO) formation. Nanoparticles were not seen in the AFM images when atomic oxygen was not present, while more particles were seen as the amount of atomic oxygen was increased. Additionally, Raman vibrational spectra indicated that oxygen was no longer bound to the graphene surface after the diethyl zinc deposition. Density functional theory (DFT) calculations were used to elucidate the mechanism by which the nanoparticles form.17 It was determined that the epoxide oxygen atom would insert between the zinc and the carbon atoms of EZE to form ethyl, ethoxy zinc (EOZE). In this process, the relatively weak zinc−carbon bond was replaced by a stronger zinc−oxygen bond, resulting in removal of the oxygen from the graphene. The energetics revealed a −3.02 eV driving force for this reaction to occur. It was also shown to be energetically favorable for EOZE to abstract a second epoxide oxygen to form diethoxy zinc (EOZOE). Further information about the initial steps in the mechanism for nanoparticle formation was determined by examining the reactions between EZE, EOZE, and EOZOE. Linear zinc−oxygen chains along with six- and eight-membered zinc−oxygen ring structures were Received: September 29, 2015 Revised: December 18, 2015
A
DOI: 10.1021/acs.jpcc.5b09511 J. Phys. Chem. C XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry C
and more complex species. We then study the formation mechanism and structure of the nanoparticles using MD simulations that are supplemented by a Monte Carlo acceleration approach for circumventing certain high energy barriers. Finally, we examine the electronic properties of the nanoparticles using time-dependent density functional theory (TDDFT) calculations.
proposed as possible products. These reactions, which expel butane and/or diethyl ether, were found to be exothermic by (approximately) 1 eV. Another DFT study examined the single and double hydrolysis of EZE as well as the dimerization and tetramerization of the resulting species.18 These DFT calculations serve to identify the initial reaction pathways and intermediates for zinc oxide nanoparticle growth. While these DFT calculations are important to understand the beginning reactions in nanoparticle formation, it is difficult to use electronic structure calculations to investigate the mechanism when the number of atoms approaches hundreds or even thousands of atoms. It is also not trivial to a priori identify all possible reactions that may play a role in nanoparticle formation. Because one would like to perform molecular dynamics (MD) simulations to sample reactions that lead to nanoparticle growth, it should be noted that it is necessary to use reactive force fields to describe bond formation and breakage. This approach has been used in past work,19−26 but force field parametrization is often complex. Also, the accuracy of reactive force fields is often disputed because they attempt to mimic complicated electronic effects using simple potential energy functions. Parameters must be determined, typically for each element and pair of elements, before a simulation can begin, and it is not always clear which parameters will be applicable to environments outside what was used in the fitting procedure. In this study, we use the reactive force field originally developed for hydrocarbons in 2001 by van Duin et al. known as ReaxFF.26 ReaxFF defines a bond between two atoms based on the “bond order”, which is a simple function of the interatomic distance. Bonds, angles, torsions, and other important interactions are then taken to be functions of this bond order. ReaxFF was first used in zinc containing systems in 2008 to model the structures and reaction dynamics for ZnO catalysts.27 The parameters were determined by fitting to Zn metal and ZnO polymorph properties as well as bond dissociation curves in Zn(OH)2. The force field was modified two years later to study water adsorption and dissociation on different ZnO crystal faces.28 ReaxFF was used later that same year to examine hydrolysis and water stability in zinc metal−organic frameworks.29 It is important to note a few things about these studies that are relevant here. In the metal−organic framework study, the authors did not include zinc−carbon parameters because these atom types were not expected to interact; thus, they deemed this interaction insignificant. In the present study it is therefore necessary to introduce zinc−carbon parameters so that EZE can be simulated. Also, in the original parametrization, the crystal properties were given a higher priority in the fitting procedure compared to the zinc−oxygen bond dissociation energy. As a consequence, the ReaxFF zinc−oxygen bond energy is underestimated by ∼20 kcal/mol. For the present work where oxygen abstraction by EZE is studied, it is important to have the correct zinc−carbon to zinc−oxygen bond energy ratio. The zinc−oxygen bond parameters in ReaxFF will therefore need to be adjusted to reflect accurately this relationship. However, this may limit the accuracy of ZnO crystal properties in the modified version of ReaxFF. The goal of this project is to develop ReaxFF to simulate ZnO nanoparticle growth from EZE, as initiated on an oxidized graphene surface. We begin by developing new ReaxFF parameters that describe the epoxide oxygen abstraction by EZE, as well as subsequent reactions involving EOZE, EOZOE,
■
METHODS Potential Energy Surface. We begin by reviewing the portions of ReaxFF relevant to this paper. Only the functions that contain parameters which were altered or added to the potential energy surface will be presented. The details of the original ReaxFF formulation can be found in ref 26. ReaxFF has been reformulated since its incarnation and we refer the reader to refs 30 and 31 for more recent explanations about the form of the potential energy surface. For consistency, the symbolic notation used in the aforementioned papers will be adopted here. The total energy is decomposed into a sum of physically relevant terms that are in turn functions of bond orders. Within ReaxFF, it is assumed that the bond order between atoms i and j is given by BO′ij = BOijσ + BOijπ + BOijππ ⎡ ⎡ ⎡ ⎛ rij ⎞ p bo,2 ⎤ ⎛ rij ⎞ p bo,4 ⎤ ⎛ rij ⎞ p bo,6 ⎤ ⎥ + exp⎢p ⎜ ππ = exp⎢pbo,1 ⎜ σ ⎟ ⎥ + exp⎢pbo,3 ⎜ π ⎟ ⎟ ⎥ bo,5 ⎢⎣ ⎥⎦ ⎝ r0 ⎠ ⎥⎦ ⎝ r0 ⎠ ⎝ r0 ⎠ ⎥⎦ ⎣⎢ ⎣⎢
(1)
where rij is the interatomic distance and the pbo’s and r0’s are fitting parameters. The prime is used to indicate an uncorrected bond order (i.e., that calculated from the simulation coordinates). It was found that summing the bond orders calculated from eq 1 would often lead to a value higher than that elements’ valency. As such, a correction with no additional parameters is applied to eq 1 to give BOij (eq 3 in ref 26). Note that as two atoms separate, their bond order smoothly goes to zero. The bonding energy between two atoms is then given by E bond = −Deσ BOijσ exp[pbe,1 (1 − (BOijσ ) p be,2 )] − Deπ BOijπ − Deππ BOijππ
(2)
where the De’s and ppe’s are fitting parameters and the BOij’s are the corrected bond orders. Because the bond orders smoothly go to zero, there is no discontinuity in the bond energy as the interatomic distances change. The valence angle potential energy is given by Eval = f7 (BOij )f7 (BOjk )f8 (Δj , pv,1 , pv,2 ){ka − ka exp[−k b(Θ0 − Θijk )2 ]}
(3)
where ka and kb are fitting parameters. Θ0 is also a fitting parameter that represents the equilibrium angle and Θijk is the angle between vectors rij and rjk. The f 7 functions ensure the valence angle energy goes to zero if either bond order goes to zero. Finally, f 8 is a function of two fitting parameters (pv,1 and pv,2) as well as the degree of over/undercoordination of the central atom j (Δj). Nonbonded van der Waals (vdW) interactions are treated between all atom pairs via the following equation ⎧ ⎡α ⎛ ⎡ ⎛ f (rij) ⎞⎤ f (rij ⎞⎤⎫ ⎪ ⎪ ⎟⎟⎥ − 2·exp⎢ vdw ⎜⎜1 − 13 ⎟⎟⎥⎬ EvdW = ϵvdW ⎨exp⎢αvdW ⎜⎜1 − 13 ⎪ ⎪ ⎥ ⎥ ⎢ r 2 r ⎢ ⎝ ⎠ ⎝ ⎠ ⎦ ⎦ ⎣ ⎣ vdW vdw ⎭ ⎩
(4) B
DOI: 10.1021/acs.jpcc.5b09511 J. Phys. Chem. C XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry C Table 1. Final Zinc−Carbon and Zinc−Oxygen Bond Order and Bond Energy Parameters Zn−C Zn−O
rσ0 (Å)
pbo,1
pbo,2
Dσe (kcal/mol)
pbe,1
pbe,2
2.0927 1.9804
−0.3467 −0.3421
4.6113 5.4933
181.3657 193.7071
−0.7035 −0.5984
1.9282 1.7527
where ϵvdW, αvdW, and rvdW are fitting parameters and f13 is a distance-dependent shielding function. ReaxFF also has terms representing torsional angles, overcoordination, undercoordination, lone electron pairs, conjugation, and Coulomb interactions. None of the parameters in these terms are altered in this study; therefore, we do not include their formulas here. As noted above, ReaxFF was previously developed for the Zn, C, O, and H system, with the exception that the zinc− carbon bonding parameters were omitted.29 For this zinc− carbon interaction, it is important to determine only the σ bond order parameters (pbo,1, pbo,2, and rσ0) and the σ bond energy parameters (Dσe , pbe,1, and pbe,2) because we do not expect zinc and carbon to form double or triple bonds. We also decided to refit the zinc−carbon van der Waals parameters because this pair can now form a σ bond. As previously stated, we will also determine new parameters for the zinc−oxygen bond energy. Valence angle parameters for C−Zn−C, Zn−C−H, Zn−C−C, C−Zn−O, and Zn−O−C were also parametrized. All other parameters are taken from ref 29. We follow as closely as possible the method by which the parameters were found previously to fit the new parameters. We first define the training set that will be used to fit the parameters by calculating various bond and angle potential energy surfaces. This consists of energies calculated using NWChem version 6.1.1 with the B3LYP functional and the 6311+G* basis set. The molecules used in the training set were EZE, EOZE, EOZOE, MZM, and MOZM where M represents methyl. Zn−C and Zn−O bond dissociation profiles were calculated by stretching the bond of interest and optimizing the remaining degrees of freedom. For each stretching potential, the zero of energy was taken to be the triplet state energy when the bond was stretched to 7 Å. The valence angle potential energy surfaces were calculated in a similar fashion where the angle of interest was fixed at a particular value while optimizing the other degrees of freedom. For each angle potential, the zero of energy was arbitrarily defined to be that of the largest angle. The algorithm used to determine the parameters is a commonly used scheme in which the parameters were iteratively optimized one at a time to best fit the training set.32 Here, we define the error function to be minimized as Error =
∑ wi(piReaxFF − piDFT )2 i
needed to perform ReaxFF simulations.33−36 Periodic boundary conditions and the NVT ensemble were used in all cases. The temperature was held at the desired temperature using the Berendsen thermostat with a coupling time of 100 fs.37 The time step of the simulations was set to 0.25 fs. In the ALD simulations, we began with a 180 carbon atom graphene sheet with two epoxide groups. This epoxide coverage is fairly consistent with that estimated in the experiments of 0.1−3% depending on the oxygen exposure. EZE was randomly placed above the graphene sheet. The system was equilibrated over 12.5 ps as the temperature was raised from 100 K to the temperature of interest. The total linear and angular momentum of the molecule in the gas phase was kept at zero to avoid deposition. After equilibration, the velocity of the zinc atom was given a random velocity based on the Boltzmann distribution of speeds and the temperature of the simulation. The sign of the velocity component normal to the surface was set to direct the zinc atom toward the epoxidized side of the graphene sheet and the ethyl groups bonded to the zinc follow. This simulation was then run for 500 ps. The output geometry could then be used as a starting point for another ALD run by adding another gas-phase molecule and more epoxide groups. Repeating this many times is one method we used to study the ZnO nanoparticle growth by ALD. As we will explain later, purely gas-phase simulations of EOZOE were also used to facilitate larger particle growth. These used the same simulation parameters as the ALD simulations with a total simulation time of 1.25 ns. Optical Spectra. The TDDFT calculations were performed using version 2013.01 of the ADF molecular modeling suite.38−42 We again used the B3LYP functional for these calculations. The TZP basis set was used for all atoms, and the frozen core approximation was used. COSMO-RS was used to incorporate the effects of a water solvent in these calculations (as is typically present in optical measurements). The error tolerance in the square of the excitation energies was set to 8 × 10−6 hartree, and 50 excitations were calculated. The peaks were broadened by a Gaussian profile with a 0.5 eV width. Excitations within 0.5 eV of the highest-energy excitation were not included in the spectral line shape.
■
RESULTS Parametrization Results. We begin by presenting the final parameters of the fitting procedure. The full parameter file used in the MD simulations is provided in the Supporting Information. The bond order and bond energy parameters for Zn−C and Zn−O pairs are given in Table 1. The Zn−O σ bond distance parameter (rσ0) was fixed at that used in ref 28. The refit Zn−O bond energy parameter (Dσe ) has been increased by more than 30 kcal/mol to improve the error in the original parametrization’s dissociation energy. The final Zn−C van der Waals parameters are shown in Table 2. The Zn−C and Zn−O bond dissociation curves are shown in Figure 1. It is perhaps unsurprising that very good agreement is seen between the ReaxFF and the DFT calculations because the weight (wi in eq 5) associated with bonds was larger than that for the valence angles. This provided control over the fitting procedure to
(5)
where the sum runs over the training set; wi is a predetermined weight corresponding to point i, and the pi’s represent the ReaxFF and DFT energy of point i (both at the DFT optimized geometry). Equation 5 is first evaluated at the current ReaxFF parameter set and then at the current parameter set with a small positive and negative change in the parameter of interest. A quadratic function is fit to these three data points, and the parameter was changed by minimizing the function inside the parameter bounds. The next parameter in the set was chosen, and the routine was repeated. This procedure was iterated until convergence was achieved or minimal changes were found when changing a parameter. MD Simulations. All simulations were performed using the LAMMPS simulation package which includes the subroutines C
DOI: 10.1021/acs.jpcc.5b09511 J. Phys. Chem. C XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry C Table 2. Final Zinc−Carbon van der Waals Interaction Parameters Zn−C
Table 3. Final Valence Angle Parameters
−1
ϵvdW (kcal/mol)
rvdW (Å)
αvdW (Å )
0.3226
1.8473
10.4746
C−Zn−C Zn−C−H Zn−C−C C−Zn−O Zn−O−C
ensure the correct Zn−C to Zn−O bond energy ratio was obtained by ReaxFF. As previously mentioned, it is important to correctly reproduce this ratio because the first proposed step in the mechanism involves the replacement of a weaker Zn−C bond with a stronger Zn−O bond. Figure S1 in the Supporting Information shows that this new ReaxFF parametrization reproduces the correct energy ordering of four ZnO crystal structures, as was previously published.27,28 The valence angle parameters are provided in Table 3. The potential energy surfaces for these fits are shown in Figures 2 and 3. ReaxFF is comparable to the DFT results with errors less than ∼4 kcal/mol in the relevant angle ranges. This is similar to the bond dissociation curves; therefore, we know the weights used are capable of obtaining a fit that balances the bonds and valence angle errors. No torsion angle parameters were fit because their potential energy surfaces (not shown) were generally featureless and the energy scales are much smaller than bonds and valence angles (
