Article pubs.acs.org/JPCC
Interaction of Alkylamines with Cu Surfaces: A Metal−Organic ManyBody Force Field Shih-Hsien Liu† and Kristen A. Fichthorn*,†,‡ †
Department of Chemical Engineering and ‡Department of Physics, The Pennsylvania State University, University Park, Pennsylvania 16802, United States S Supporting Information *
ABSTRACT: Hexadecylamine (HDA) and alkylamines, in general, are key molecules in the shape-selective synthesis of Cu nanostructures. To resolve certain aspects of these syntheses, we develop a classical many-body force field to describe the interactions of HDA with Cu surfaces. We parametrize the force field through force and energy matching to results from first-principles density functional theory (DFT). Our force field reproduces the DFT binding energies and configurations of self-assembled HDA layers on Cu(100) and Cu(111) at various coverages. We implemented the force field in classical molecular dynamics (MD) simulations to resolve various HDA self-assembled-layer structures on Cu(100) in vacuum, and we find that HDA layers undergo a continuous structural transition through various ordered layers at high coverage to disordered layers at lower coverages. We probed pentylamine (PA), decylamine (DA), and HDA binding on Cu surfaces in vacuum with MD and find that DA forms self-assembled layers, but PA layers disorder at experimental temperatures. We investigated HDA monolayers on Cu surfaces in an aqueous medium with MD and found that the self-assembled monolayer structure in vacuum is retained. The long and hydrophobic alkyl tails in the self-assembled HDA monolayer repel water molecules and would prevent Cu oxidation, which agrees with experiment.
■
INTRODUCTION The synthesis of Cu nanostructures is of considerable current interest for applications such as flexible nanowire electrodes in transparent conducting films,1−13 conductive ink,14 and catalysis.15−20 In these applications, the unique properties of these nanocrystals are significantly influenced by their size and shape. To this end, alkylamines have been widely used as capping agents for the solution-phase synthesis of various Cu nanostructures. These alkylamines include octylamine,21,22 dodecylamine,21 tetradecylamine,12 octadecylamine,9,15,23,24 oleylamine (OLA), 2,4,11,16,23,25,26 and hexadecylamine (HDA).1,3,6,13−15,19,20,22,23,27−31 The role of these molecules in achieving controlled nanocrystal shapes may be multifold: They deter oxidation of the Cu nanostructures,1,6,13,28,30 prevent random aggregation, and direct the shapes that form by influencing either the kinetics or the thermodynamics of these structures. A detailed understanding of the function of these molecules is important in achieving controlled and reproducible products. Atomic-scale theory and simulations can help to propel our understanding of the solution-phase, shape-selective growth, and unique properties of nanostructures. First-principles studies can resolve interactions that are important in describing the growth and properties of metal nanocrystals.32−39 However, these studies are limited to relatively small systems and are typically conducted in a zero-temperature vacuum environment that is far removed from the solution environment of nanostructure synthesis. On the contrary, classical molecular © XXXX American Chemical Society
dynamics (MD) simulations allow for faithful recreation of the solution-phase environment, in many cases on length and time scales relevant to understanding nanocrystal growth, provided that these can be based on reliable interatomic force fields. Using MD, progress can be made in understanding the solution-phase binding of capping agents to nanocrystal surfaces,40−42 solid−liquid interfacial free energies and nanocrystal Wulff shapes,43 as well as kinetic structures and the kinetics of nanocrystal growth.44−48 In this work, we enable such studies by presenting a force field designed to describe the binding of alkylamines to Cu surfaces. An important aspect of a force field for predicting shapeselective nanocrystal synthesis is the capability to reproduce facet-selective binding consistent with experimental studies. For example, recent experimental studies show that HDA1,3,15,19,20,27,29,30 and OLA10,16,25,26 tend to promote the formation of {100}-faceted Cu nanostructures, such as nanowires and nanocubes, in solution-phase syntheses. These results suggest that HDA and OLA bind selectively to the {100} facet, as the {111} facet has the lowest energy for bare Cu,39 and we would expect primarily {111}-faceted nanostructures to be preferred thermodynamically.43 Our recent studies with dispersion-corrected density functional theory (DFT) also indicate that HDA binds preferentially to Cu(100) surfaces.39 In these studies, we found that HDA forms self-assembled Received: August 7, 2017 Published: September 15, 2017 A
DOI: 10.1021/acs.jpcc.7b07861 J. Phys. Chem. C XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry C
Figure 1. Top-down and side views of the DFT-optimized HDA structure in the unit cell and underlying Cu surface atoms for (5 × 3)-Cu(100) with (a) uniform γ and (b) mixed γ, for (c) c(2 × 6)-Cu(100), as well as for (d) p(2 × 2)-Cu(100).39 The red rectangle is the unit cell of the pattern. Nitrogen is blue, carbon is black, and hydrogen is white.
monolayers on Cu(100) and Cu(111), similar to those observed experimentally for alkanethiols on Cu(111)49 and Au(111).50−53 Interestingly, we found that the binding energy of HDA molecules in the (5 × 3) overlayer identified for Cu(100) was higher than that of molecules in the ( 3 × 3 )R30° overlayer that we found for Cu(111).39 To develop a force field consistent with DFT and experiment, we require preferential binding of HDA molecules to the Cu(100) surface, which has a more open structure than Cu(111). A similar scenario occurs for polyvinylpyrrolidone (PVP), which was found in first-principles DFT calculations to bind more strongly to Ag(100) than to Ag(111)33,36 and to Pt(100) than Pt(111).38 It is difficult to reproduce the {100} binding preference with pair potentials54 and we developed a Metal−Organic Many-Body (MOMB) force field46,55 to model this preference for the PVP/Ag system. In the alkylamine/Cu system, there is charge transfer (in both DFT39 and experiment9) from the amine group of HDA to Cu surface atoms, which plays an important role in dictating {100} facet selectivity; a similar scenario occurs for the O atom in the PVP/ Ag system.33,36 As we elaborated elsewhere46 and we will discuss below, the capability of the MOMB force field to mimic this charge transfer and its ramifications for metal−organic binding conformations is important to its success.
Figure 2. Top-down and side views of the DFT-optimized HDA structure in the unit cell and underlying Cu surface atoms for ( 3 × 3 )R30°-Cu(111) with (a) uniform γ and (b) mixed γ as well as for (c) p(2 × 2)-Cu(111).39 The red rectangle is the unit cell of the pattern. Nitrogen is blue, carbon is black, and hydrogen is white.
■
Table 1. Morse Potential Parameters D0, α, and r0 in Equation 1 and the Scaling Factors for One-Way Amine H → Cu ( f H) and N → Cu (f N) Electron-Density Functions in Equations 2 and 3
METHODS DFT Calculations. We obtained force-field parameters by fitting to results of HDA binding on Cu surfaces from firstprinciples DFT calculations. We use the Vienna Ab Initio Simulation Package (VASP)56−58 employing projector augmented-wave (PAW) potentials59 to describe electron−ion core interactions and the generalized-gradient approximation (GGA) exchange-correlation functional by Perdew, Burke, and Ernzerhof (PBE)60 for our DFT calculations. Long-range van der Waals (vdW) interactions are included using the DFT-D2 method of Grimme,61 where we implement a dispersion coefficient (C6) and vdW radius (R0) derived by Ruiz et al. for B
interaction
D0 (eV)
α (Å−1)
r0 (Å)
f H, f N
amine H−Cu aliphatic H−Cu C−Cu N−Cu
0.00176 0.00157 0.00103 0.00418
1.87 1.64 1.35 4.13
2.51 2.82 4.27 2.31
1.29
4.24
DOI: 10.1021/acs.jpcc.7b07861 J. Phys. Chem. C XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry C
predicted by the EAM potential (3.615 Å65,66) agrees well with experiment,67 as do the cohesive energy, bulk modulus, and shear modulus.65 To model the HDA intra- and intermolecular interactions, we use the CHARMM General Force Field (CGenFF)68,69 implemented in the CGenFF program70 that performs atom typing and assignment of bond, angle, and dihedral parameters as well as charges by analogy.71,72 We also use the CHARMM36 force-field files73 to assign Lennard-Jones parameters to each type of atom in an HDA molecule. We note that the CGenFF program74 and the CHARMM36 general force field75 have been applied to small organic molecules adsorbed on carbon nanotubes in MD simulations. We use the framework of the MOMB force field46,55 to model HDA−Cu interactions. In this potential, atoms (A) in HDA interact with Cu surfaces mainly via a pair potential (ϕA−Cu) defined as
Table 2. HDA Binding Patterns on Cu(100) and Cu(111) from DFT39 and Their Corresponding Properties Defined in the Text pattern
ρHDA (Å−2)
G
x axis
y axis
0.061 0.051 0.038 0.059
6 4 1 2
[011] [011] [011] [2̅11]
[01̅1] [01̅1] [01̅1] [01̅1]
0.044
2
[1̅10]
[1̅1̅2]
(5 × 3)-Cu(100) c(2 × 6)-Cu(100) p(2 × 2)-Cu(100) ( 3 × 3 )R30°-Cu(111) p(2 × 2)-Cu(111)
ϕA−Cu(rij) = D0,A−Cu[e−2αA−Cu(rij − r0,A−Cu) − 2e−αA−Cu(rij − r0,A−Cu)] − s6fdamp (rij , R 0,A , R 0,Cu)C6,A−Curij−6
(1)
where rij is the distance between atomic species i and j. The first term in eq 1 is a Morse potential with parameters D0,A−Cu, αA−Cu, and r0,A−Cu to account for the short-range and directbonding interactions. The second term in eq 1 accounts for long-range vdW interactions in the form prescribed by Grimme,61 which has been implemented in the USER-MISC package of LAMMPS.46,55,76 Here s6 is a global scaling factor, fdamp is a damping function associated with rij, the vdW radii for atoms A and Cu are R0,A and R0,Cu, and C6,A−Cu is the dispersion coefficient for the A−Cu interactions. The values of C6 and R0 for atomic species in HDA were given by Grimme,61 and those for Cu were given by Ruiz62 to account for bulk screening effects. We use a cutoff radius of 6.0 and 12 Å for the Morse potential and vdW interactions, respectively. As discussed above, a pair potential does not adequately capture the surface-selective binding observed for the HDA/Cu system. We know from prior studies46,55 that we can better match DFT results by including a term in the molecule−surface potential to interface with the EAM potential. This one-way electron density function mimics charge transfer from the amine group in HDA to the Cu surface, consistent with our findings in a recent DFT study of this system.39 In the EAM potential for Cu, we include the electron densities supplied one-way from N to Cu (ρN→Cu) and from amine H to Cu (ρH→Cu). The functions are defined as
Figure 3. HDA molecule in Cartesian coordinates with angles θ, β, and γ, as defined in the text. The x−y plane is the surface plane, and nitrogen is placed at the origin as the reference point. Nitrogen is blue, carbon is black, and hydrogen is white.
Cu to account for bulk screening effects.62 Details of the DFT calculations are included in the Supporting Information. Details of the Force Field. We model the binding of HDA on Cu surfaces within the LAMMPS MD simulation code.63,64 To model interatomic interactions, we use three different types of potentials. We use an embedded-atom method (EAM) potential65,66 for Cu−Cu interactions. The Cu lattice constant
ρN → Cu (rij) = fN ρCu−Cu (rij)
(2)
ρH → Cu (rij) = fH ρCu−Cu (rij)
(3)
Table 3. Force-Field and DFT39 (in parentheses) Properties Defined in the Text for the Three HDA Binding Patterns on Cu(100)a γ uniform mixed
a
Ebind (eV)
error
1.96 (1.93) 1.87 (2.01)
1.44% 6.60%
1.94 (1.96)
1.38%
1.93 (1.96)
1.47%
EHDA−HDA (eV) (5 × 3)-Cu(100) 1.46 (1.46) 1.37 (1.52) c(2 × 6)-Cu(100) 1.31 (1.41) p(2 × 2)-Cu(100) 1.07 (1.24)
EHDA−Cu (eV)
⟨dN−Cu⟩ (Å)
⟨θ⟩
0.50 (0.48) 0.50 (0.49)
2.33 (2.25) 2.33 (2.27)
0.8° (1.7°) 0.8° (1.8°)
0.63 (0.55)
2.20 (2.23)
31.7° (35.9°)
0.86 (0.72)
2.14 (2.25)
50.6° (50.8°)
Corresponding DFT-optimized structures are shown in Figure 1. C
DOI: 10.1021/acs.jpcc.7b07861 J. Phys. Chem. C XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry C
Table 4. Force-Field and DFT39 (in parentheses) Properties Defined in the Text for the Two HDA Binding Patterns on Cu(111)a γ
Ebind (eV)
uniform mixed
a
EHDA−HDA (eV)
error
1.86 (1.83) 1.84 (1.87)
1.45% 1.68%
1.82 (1.89)
3.86%
EHDA−Cu (eV)
(√3 × √3)R30°-Cu(111) 1.43 (1.39) 0.42 (0.44) 1.41 (1.45) 0.43 (0.42) p(2 × 2)-Cu(111) 1.20 (1.28) 0.62 (0.62)
⟨dN−Cu⟩ (Å)
⟨θ⟩
2.40 (2.30) 2.41 (2.28)
2.6° (2.6°) 0.3° (0.2°)
2.33 (2.25)
38.5° (38.8°)
Corresponding DFT-optimized structures are shown in Figure 2.
Figure 4. Comparison of total binding energies per HDA molecule predicted by the force field (FF) to those by DFT for the five HDA binding patterns in Table 239 and variants of these patterns given in the Supporting Information. The solid line indicates a prefect fit of the FF to DFT and the dotted lines indicate a 10% error.
Figure 5. Radial distribution functions (g(r)) for nitrogen atoms projected onto the x−y plane with standard error bars for the three HDA binding patterns on Cu(100) with high HDA packing density and the (19 × 3)-Cu(100) shown in Table 5 at a temperature of 373 K. Vertical lines show g(r) for (5 × 3)-Cu(100) at 0 K as a reference.
where the scaling factors f N and f H are used to adjust the amount of electron density supplied from N to Cu and from amine H to Cu, respectively, and ρCu−Cu is the electron density between Cu atoms. Thus a modified EAM potential for Cu (ϕCu) is given as
atoms. The superscript on the sum indicates the sum runs over the specified species or species pairs. We note that FCu, ρCu−Cu, and ϕCu−Cu were given by the original EAM potential for Cu,65,66 and the modified EAM potential was implemented in LAMMPS via the Finnis−Sinclair option.77 Force-Field Fitting. For HDA−Cu interactions described in the pair potential in eq 1, we note that D0,A−Cu, αA−Cu, and r0,A−Cu are unknown parameters, and s6 may need to be adjusted to better estimate vdW interactions.37,55 We determined these parameters for each of four designated atom types (A) in an HDA molecule: amine H bonded to N, aliphatic H bonded to C, aliphatic C, and N. Additionally, to mimic the effect of charge transfer from the amine group to the surface, we determined the values of f N and f H for the amine H and the N in eqs 2 and 3.
Cu
ϕCu =
i
H − Cu
+
⎛ Cu − Cu
∑ FCu⎜⎜ ∑ ∑ i,j
⎝
N − Cu
ρCu−Cu (rij) +
i≠j
⎞ 1 ρH → Cu (rij)⎟⎟ + 2 ⎠
∑
ρN → Cu (rij)
i,j Cu − Cu
∑
ϕCu−Cu(rij)
i≠j
(4)
where FCu is the Cu embedding energy, which is a function of electron density ρ and ϕCu−Cu is the pair potential between Cu
Table 5. Proposed HDA Binding Patterns on Cu(100) and Cu(111) in This Study and Their Corresponding Properties Defined in the Text pattern (13 × 3)-Cu(100) (5 × 3)-Cu(100) (12 × 3)-Cu(100) (19 × 3)-Cu(100) (26 × 3)-Cu(100) (33 × 3)-Cu(100) (7 × 3)-Cu(100) (16 × 3)-Cu(100) (17 × 3)-Cu(100) ( 3 ×
3 )R30°-Cu(111)
ρHDA (Å−2) 0.063 0.061 0.060 0.059 0.059 0.059 0.058 0.057 0.054 0.059
(high) (high) (high) (medium) (medium) (medium) (low) (low) (low) (medium) D
G
x axis
y axis
Ebind,MD (eV)
16 6 14 22 30 38 8 18 18 2
[011] [011] [011] [011] [011] [011] [011] [011] [011] [2̅11]
[011̅ ] [01̅1] [011̅ ] [01̅1] [01̅1] [01̅1] [011̅ ] [01̅1] [01̅1] [01̅1]
1.84 1.87 1.88 1.88 1.88 1.88 1.88 1.87 1.75 1.81
± ± ± ± ± ± ± ± ± ±
0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01
DOI: 10.1021/acs.jpcc.7b07861 J. Phys. Chem. C XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry C
(15 × 12 × 26) and (7 3 × 12 × 23) supercells for Cu(100) and Cu(111), respectively, where the unit to describe the x and y supercell axes in the surface plane is the nearest-neighbor distance and that in z direction normal to the surface is the interlayer spacing. The HDA/Cu system has six Cu layers, where the bottom three layers were fixed at the bulk termination and all of the other atomic coordinates were fixed at DFT-optimized configuration. From these configurations, we generated 15 different HDA conformations on each Cu surface by translating the molecules along the z axis normal to the Cu surface. We translated the molecules up to 0.9 Å above and down to 0.5 Å below the original conformation, in 0.1 Å increments to probe both attractive and repulsive forces exerted by the Cu surface on each atomic species in the molecules. For each fixed HDA conformation, both DFT and force-field calculations were performed to obtain the force and energy profiles for HDA−Cu interactions. The force profiles were obtained using
Figure 6. Similar to Figure 5, except for medium HDA packing density shown in Table 5.
HDA−Cu total HDA Fmlk = Fmlk − Fmlk
(5)
HDA−Cu Fmlk
where is the k Cartesian component of the force exerted by the Cu surface on atom l in the HDA conformation HDA total is the total force, and Fmlk is the intra- and m, Fmlk intermolecular force, which was calculated using the same HDA conformation in a supercell with the same size as for the HDA/ Cu system, only without Cu. The energy profiles were obtained using EmHDA−Cu = EmHDA + EmCu − Emtotal
(6)
where EmHDA−Cu is the interaction energy between the HDA conformation m and the Cu surface, EmHDA is the energy of the same HDA conformation in a supercell with the same size without Cu, EmCu is the energy of the Cu surface in a supercell with the same size without HDA, and Emtotal is the energy of the HDA/Cu system. The initial values of unknown variables (D0,A−Cu, αA−Cu, r0,A−Cu, f N, f H, and s6) were chosen according to our previous study,55 and they were adjusted randomly using a Python-
Figure 7. Similar to Figure 5, except for (33 × 3)-Cu(100) and low HDA packing density shown in Table 5.
We determined the force-field parameters using the most favored HDA binding patterns on Cu(100) and Cu(111) (cf., Figures 1a and 2, respectively) in our recent DFT study.39 To carry out the parametrization, we replicated these structures in
Figure 8. Top-down and side views of the equilibrated (a) PA, (b) DA, and (c) HDA configurations for (5 × 3)-Cu(100) at 373 K. Nitrogen is blue, carbon is black, and hydrogen is white. E
DOI: 10.1021/acs.jpcc.7b07861 J. Phys. Chem. C XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry C
Figure 9. Top-down and side views of the equilibrated (a) PA, (b) DA, and (c) HDA configurations for ( 3 × Nitrogen is blue, carbon is black, and hydrogen is white.
Here, the first, second, and third terms correspond to the force, energy, and energy difference between the two Cu surfaces, respectively. m = 1 corresponds to the HDA conformation closest to the surface and the conformation furthest from the surface is M = 15. L = 24 is the number of atoms per HDA molecule within the vdW cutoff to the surface atoms, and the weighting factor w = L−1 accounts for the fact that the number of energy observations is less than that of forces. In the algorithm, we generate random trial parameter values, and only the values that lead to a cost function smaller than the current value are accepted for the next iteration. When the value of the cost function converged to a minimum, s6 = 0.60, and the other variables are summarized in Table 1. Force-Field Testing. Following parametrization of the force field, we tested its accuracy by performing energy minimization on various DFT-optimized HDA structures on Cu(100) and Cu(111). For an HDA molecule in the gas phase, we used a cubic supercell with a side length of 50 Å. The two Cu surfaces, Cu(100) and Cu(111), with adsorbed HDA were modeled using supercells with six Cu layers, where the bottom three layers were fixed at the bulk termination and all of the other atomic coordinates were relaxed. For the HDA/Cu system in vacuum, we used a supercell with length of 70 Å in the surface normal. We performed energy minimization using the FIRE algorithm78 until the maximum force on any Cartesian component of any relaxed atom (i.e., f max in LAMMPS) was
