ARTICLE pubs.acs.org/Langmuir
Determination of Binding Energy and Solubility Parameters for Functionalized Gold Nanoparticles by Molecular Dynamics Simulation Brian J. Henz,* Peter W. Chung, Jan W. Andzelm, Tanya L. Chantawansri, Joseph L. Lenhart, and Frederick L. Beyer U.S. Army Research Laboratory, Aberdeen Proving Ground, Maryland 21005, United States ABSTRACT: The binding energy, density, and solubility of functionalized gold nanoparticles in a vacuum are computed using molecular dynamics simulations. Numerous parameters including surface coverage fraction, functional group (CH3, OH, NH2), and nanoparticle orientation are considered. The analysis includes computation of minimum interparticle binding distances and energies and an analysis of mechanisms that may contribute to changes in system potential energy. A number of interesting trends and results are observed, such as increasing binding distance with higher terminal group electronegativity and a minimum particleparticle binding energy (solubility parameter) based upon surface coverage. These results provide a fundamental understanding of ligand-coated nanoparticle interactions required for the design and processing of highdensity polymer composites. The computational model and results are presented as support for these conclusions.
’ INTRODUCTION There is currently significant interest in the processing and manufacture of polymer composite materials with nanometersized fillers, due to the potential increases in mechanical, chemical, and electrical properties they can provide.1,2 The addition of nanoparticles (NPs) in polymer matrices has been shown to drastically alter the properties of the composite material over composites produced with micrometer sized particles.3,4 One advantage of nanometer-sized fillers, or in this effort NPs, is that they have a large surface area to volume ratio. The high surface area to volume ratio of nanometer-sized fillers decreases the amount of filler material required to produce a property change in the composite while minimizing the possibly deleterious effect on the desirable properties of the polymer matrix such as flexibility. Experimental observations have shown that controlled, regularly spaced, and dense NP packing can lead to better control of polymer composite properties of interest and offer improved results over randomly distributed loadings.5 One method used to control the placement and distribution of NPs is to use block copolymers as the polymer matrix as opposed to a homogeneous material.3,68 The efficient use of block copolymers requires a fundamental understanding of NPNP and NPmatrix interactions. Multiscale simulation methods are used to investigate phenomena that occur on multiple time and length scales. Nanoparticle composite materials, by their nature, fit well into this framework. The multiscale modeling method used for this work couples the atomic and meso length scales through the development of an empirical potential for NPNP and NPmatrix interactions. This coupling can be performed in multiple ways r 2011 American Chemical Society
from simple parameter development to tightly coupled simulations that interact iteratively. Here we have started by computing a single parameter, namely, the solubility parameter,9 of functionalized NPs to describe NPNP and NPpolymer interactions. The solubility parameter is dependent upon the cohesive energy of the pure materials and takes into account binding energy and packing density or bond length. The work detailed here is an investigation of the surface interactions between NPs and the effects of the ligand coating on those interactions. In this paper, we compute the interactions of functionalized NPs in a vacuum using molecular dynamics simulations. Analysis of these simulations provides qualitative results and insight into the interactions of functionalized NPs. We have achieved an understanding of the NPNP interactions of functionalized NPs through the systematic modeling effort described here. The functionalized NPs considered in this effort are alkanethiol ligand-coated gold NPs. Gold NPs are used as a model system here, as they possess interesting properties including electrical, magnetic, optical, and physical properties.10 In addition, there exists a wealth of experimental11,12 and computational data13,14 including simulations of NPs with explicit solvents.15 These alkanethiol ligands often form a self-assembled monolayer (SAM) on the gold surface. SAMs composed of polymer ligands are used to passivate NP surfaces, lower friction coefficients, and modify the chemical properties of NPs.10,16 In this effort, we have focused on the controlled distribution of functionalized NPs Received: February 7, 2011 Revised: April 13, 2011 Published: May 17, 2011 7836
dx.doi.org/10.1021/la2005024 | Langmuir 2011, 27, 7836–7842
Langmuir
ARTICLE
Table 2. Nonbonded Interaction Potential Parameters OPLS Parameters εLJ (kcal/mol)
σ (Å)
S
0.250
3.550
0.000
CH2
0.118
3.905
0.000
CH3
0.207
3.775
0.000
O
0.170
3.070
0.700
HOH
0.000
0.000
0.435
CH3(OH)
0.207
3.775
0.265
N
0.170
3.250
0.850
HNH2 Au
0.000 0.039
0.000 2.935
0.425 0.000
interaction site
Figure 1. Schematic of the alkanethiolate ligand, S(CH2)8X, used in this work, where X represents the ligand functional group.
Table 1. Bonded Interaction Potential Parametersa
q (e)
bond stretch Kr (kcal/mol Å2)
123 parameters
r0 (Å) pair
C12 10
-12
(Å kcal/mol) C3 10-3 (Å12 kcal/mol) r0 (Å) 12
SC CC
821 700
1.820 1.530
AuCH2
55 645
33.98
0.86
CN
490
1.335
AuCH3
67 768
41.34
0.86
CO
317
1.522
OH
317
1.000
NH
434
1.010
Morse potential
angle Kθ (kcal/mol rad2) SCC
203 876
θ0 (deg) 114.4
CCX
203 876
109.5
CNH HNH
SHAKE SHAKE
119.8 120.0
COH
SHAKE
108.5
torsion A1
(kcal/mol) 2.218
A2
2.905
A3
3.136
A4
0.731
A5
6.272
The angle parameters with SHAKE are fixed at the specified angle using the SHAKE method.24 a
using triblock copolymers. The solubility of the NPs is controlled by modifying the chemistry of the alkanethiol functional group. By controlling the solubility of the NPs, it is anticipated that the distribution of NPs will closely follow the morphology of the triblock copolymer matrix. The primary computed property of interest here is the solubility of the NPs as compared to the solubility of the phases of the triblock copolymer. By computing the solubility of functionalized NPs, it will be possible to predict their compatibility with the phases of a block copolymer.
’ SIMULATION DETAILS The molecular dynamics method as implemented in LAMMPS17 (Large-scale Atomic/Molecular Massively Parallel Simulator) is used in this work to study the dynamic properties of NP interactions. The ligands considered in this effort consist of a C8 alkylthiol chain with CH3, OH, or NH2 terminating or functional groups, that is,
pair
D0 (kcal/mol)
R (Å-1)
r0 (Å)
AuS
9.04
1.378
2.903
S(CH2)8X (Figure 1). A previous comprehensive analysis that focused on CH3 functionalized Au NPs details the EAM (embedded atom method) interatomic potential13,14 used for the NP gold core, and the bonded and nonbonded potentials for the CH3 terminated alkanethiol ligand.18 This previous effort focused on single particle dynamics with analysis of NP stability, ligand length, temperature effects, and so forth. From this analysis, the author determined the necessity of simulating the entire NP system including the Au core using a dynamic model. This is in contrast to previous efforts that assumed the Au substrate to be a flat surface1921 or a nanocrystallite.10,22 The empirical potentials here include the valence terms whose parameters are given in Table 1 and include bond stretching, bending, and torsion potentials. The united atom (UA) method23 is used to model the ligands with the exception of the terminal groups, which are modeled as fully atomistic. The functional groups that require charge considerations, namely, OH and NH2,19 are modeling fully atomistically. Bond stretching is modeled by a harmonic potential, eq 1. Ebond ¼ Kðr r0 Þ2
ð1Þ
Bending and torsion are also described by harmonic potentials as given by eqs 2 and 3, respectively. Eangle ¼ Kðθ θ0 Þ2 Etorsion ¼
∑
n ¼ 1:5
An cosn 1 ðjÞ
ð2Þ ð3Þ
The torsion parameters given in Table 1 are for the alkanethiol ligand, XCH2CH2X, where X can represent any of the united atom groups, that is, S, CH3, N, or O. The nonbonded ligand interactions are modeled using the UA OPLS (optimized potential for liquid simulations)2527 and a 123 potential (Au-CH2 and Au-CH3 only), eqs 4 and 5, respectively, with the parameters given in Table 2. The OPLS potential is modeled 7837
dx.doi.org/10.1021/la2005024 |Langmuir 2011, 27, 7836–7842
Langmuir
ARTICLE
as a combination of a van der Waals LJ (LennardJones) potential term and an electrostatic (Coulombic interactions) term, eq 4. " 6 # qi qj σ 12 σ ð4Þ þ EOPLS ¼ 4εLJ r r εr
E123 ¼
C12 C3 12 ðr r0 Þ ðr r0 Þ3
h i EMorse ¼ D0 e2Rðr r0 Þ 2eRðr r0 Þ
ð5Þ
ð6Þ
In eq 4, εLJ refers to the LJ energy parameter and ε is the dielectric constant. The LJ parameters for S, Au, and CH2 are given for reference, as they are used to compute the cross interaction parameters using the LorentzBerthelot rule. The gold substrate consists of an equilibrated NP with a diameter of 5 nm and 4093 gold atoms. Equilibration of the substrate is achieved by repeatedly heating the bare gold NP to above its size dependent melting point around 1000 K18 and slowly cooling over 1 ns by velocity rescaling back to 300 K. This results in a faceted gold core with large flat faces.28 Next, alkanethiol chains are randomly distributed over the gold surface at densities of up to 15.4 Å2 per ligand.29 After adsorption, the ligands are heated along with the gold core from 100 K up to 300 K and thermalized for 2 ns, allowing for diffusion over the gold surface so that the chains bind at suitable locations. In this work, thermalization is defined to occur when the time averaged system potential energy is constant with no applied thermostat. Finally, a constant energy (NVE) simulation is performed to check for energy conservation and a constant average temperature. The final result is a model of a functionalized gold NP, similar to those produced experimentally in the laboratory. The next section will detail how these NPs are used in simulating the NPNP interactions. One Particle Computational Procedure. As a baseline for simulating the NPNP interactions, we initially established the reference case of a single particle in vacuum. This single NP model is subsequently used to populate a larger model containing multiple NPs. The single NP simulations are parametrized by NP surface coverage and ligand functional group. Four surface coverage fractions, namely, 20%, 50%, 80%, and 100%, were chosen for analysis. In addition to the surface coverage fractions considered, three functional groups were selected for investigation, namely, CH3, NH2, and OH. In ref 18, the authors performed a comprehensive investigation of a single functionalized gold NP in vacuum. This investigation included the pair correlation function, diffusion, radial pressure and density distributions, and phase behavior for ligands capped with a CH3 functional group.18 The additional functional groups considered here, namely, OH and NH2, are not rigorously analyzed as single NPs, but many of the previous results are expected to remain qualitatively consistent. This is in spite of the fact that these polar functional groups are modeled using the all-atom method. The additional functional end groups of OH and NH2 are chosen for this effort because they represent polar hydrophilic and hydrophobic end groups, respectively, and they are common options for synthesis of functionalized NPs. For each functional group, the considered surface coverage fractions are 20%, 50%, 80%, and 100%, resulting in 12 unique parametrizations. The observed phases for these surface coverage fractions fall within the four distinct phases identified elsewhere.18,30 These phases are a striped phase at low surface coverage, where the chains lie flat along the gold surface, and intermediate structures for higher surface coverage, where some chains are partially standing and others are lying flat. At high surface coverage, the alkanethiolate chains are in either a highly structured c(4 2) phase30 at low temperatures or a more random liquid phase at higher temperatures.
Figure 2. Thermalized ligand-coated gold NPs with hydrophilic (a) OH end group, hydrophobic (b) CH3, and hydrophobic (c) NH2. The gold core is dark gray, with yellow sulfur head groups and purple CH2 backbones. For the functional groups, oxygen atoms are red, nitrogen is blue, hydrogen is white and CH3 is green. From previous efforts,30 it is expected that the hydrophilic alkanethiol ligand will extend normal to the NP surface whereas the hydrophobic ligand is more likely to bend back down to the Au surface. This bonding has the effect of decreasing the ligand coating thickness with hydrophobic end groups. The thermalized ligand coated NPs with 100% surface coverage are shown in Figure 2. A bent alkanethiol ligand is observed with the hydrophobic NH2 functional group, whereas the hydrophilic functionalized ligand with OH end groups contains much straighter chains. The hydrophobic CH3 functionalized ligands do not appear as bent as the NH2 ligands but from the core separation results presented later it is inferred that the CH3 functionalized ligands are between the OH and NH2 functionalized ligands in reference to their structure. Two Particle Computational Procedure. For each surface coverage fraction and functional group considered, six random NP orientations were generated from the thermalized models. Multiple orientations are used for consideration of multiple contact faces and to provide some variability in the computed results. The average results for these trial systems are reported in the following analysis. Here we have focused on predicting trends rather than full statistical results, which would require a much larger data set for analysis. Future efforts will duplicate these simulations many more times with different configurations in order to obtain a lower sensitivity to NP orientation. The total initial NPligand material system potential energy in the absence of a solvent is computed by summing the potential energy of two randomly chosen thermalized NP orientations. Velocity rescaling is used to maintain a constant temperature during these NVT simulations of NP approach and separation. With all orientations originating from the same thermalized model, the initial energy is constant for a specific set of parameters. After initialization, the NPs are placed together in a nonperiodic simulation box, and the nonperiodic boundary conditions are used to minimize any external long-range interactions. The initial distance between the NPs is increased until the NPNP interaction is less than 0.1% of the sum of the potential energy of two thermalized NPs. One of the NPs is then fixed by immobilizing a 2 nm diameter core of atoms in the NP. The other NP is given a small velocity, about 10% 7838
dx.doi.org/10.1021/la2005024 |Langmuir 2011, 27, 7836–7842
Langmuir
ARTICLE
Figure 3. Plot showing typical total system potential energy versus NPNP separation with 20% (a), 50% (b), 80% (c), and 100% (d) surface coverage. The arrows have been added to indicate the approach phase, corecore separation decreasing, and the separation phase. An equilibrium separation and system energy are reached at between 4.7 and 6.5 nm depending upon surface coverage, at which point the NPs are pulled apart. of the thermal velocity, directed toward the first, fixed, NP. As the NPs begin to interact through long- and then short-range forces, the system potential energy decreases and reaches a minimum at the equilibrium separation distance. After the relative motion of the NPs is determined to be zero, a thermalization step is performed prior to NP separation. Following thermalization, the NPs are separated by selecting a 2 nm diameter core of the second, mobile, nanoparticle and moving it slowly away from the fixed NP. This separation process is accomplished in 0.025 nm steps with 10 ps thermalization increments between displacements.
’ RESULTS AND DISCUSSION Many results have been computed in this effort and are discussed below including analysis of the NPNP binding energy, the Hildebrand solubility parameter, and the development of empirical potential parameters for multiscale coupling efforts. Each of these results is detailed in the following sections along with a discussion of the observations and analysis of the results. Two Particle Simulation Results. In Figure 3, the system potential energy is plotted for each of the surface coverage fractions with the NH2 functional group ligand-coated NPs. Although similar data has been collected for the CH3 and OH functional groups, it is not all shown here in the interest of brevity. The solid curves in Figure 3 show the potential energy during the approach phase, or when the NPs are moving toward each other. The separation phase curves plot the system potential energy while the NPs are being pulled apart. During the approach phase, the two NPs initially impact compressing the ligand coating and then spring back, as observed by the back and forth data line near the minimum corecore distance. This is particularly pronounced in the high, 80100%, surface coverage
models. As the ligands move and become entangled, the NPs continue to slowly approach each other while the oscillations from the initial impact approach zero. The particle separation distance eventually reaches a minimum and the average approach velocity becomes zero. Once this velocity becomes zero, the equilibrium separation distance is computed based on the NP centers of mass and the minimum potential energy is also recorded. The NPs are then separated using the method described above to calculate the unrecovered part of the system potential energy. We attribute this change of the system potential energy to path-dependency of particle interactions, which forms the basis of the unrecovered system energy. The minimum system potential energy is then subtracted from the initial, separated, system potential energy in order to compute the binding energy between the two ligand-coated NPs. The binding energy and minimum core distance are averaged over three trial simulations. The average minimum core distance and energy results are plotted in Figure 4. Changes in System Potential Energy. In Figure 4, there are two energy curves, namely, the total system potential energy computed while the NPs are approaching each other and the other computed while the NPs are being separated. The curve representing the system energy during the approach phase of the simulations is higher than the system potential energy curve computed during separation. This is an interesting trend observed in Figure 4 for all functional groups considered, as it identifies a change in system potential energy. This change is observed consistently during each simulation. For low surface coverage, that is, 20%, the approach and separation curves are similar in shape with a fairly small change in system energy. As the amount of NP surface coverage increases, the slope changes 7839
dx.doi.org/10.1021/la2005024 |Langmuir 2011, 27, 7836–7842
Langmuir
ARTICLE
Figure 4. Plots showing minimum corecore separation (a) and binding energy (b) versus surface coverage and ligand end group.
dramatically and a large drop in potential energy is observed. The likely cause of this is therefore the ligands on the NP surface. Prior to the NPNP approach simulations, the NPs are thermalized and the system potential energy is recorded. After the NPs approach each other and upon contact, the ligands form a passivating layer between the NPs and slow forward movement; this decreases the approach velocity significantly. While forward movement of the mobile NP is halted, the system potential energy continues to decrease. This occurs because the organization of the ligands is changing to accommodate the new interactions. After thermalization this reorganization stalls and the minimum NPNP separation and system potential energy are recorded. When the reorganized ligands are separated during NPNP separation, they do not reverse the reorganization process and stay at a lower energy configuration. After separation to the starting distance, the system potential energy does not return to its higher initial value and we are left with an amount of unrecovered potential energy. The binding energy is extracted into AuAu, ligandligand, and Auligand interactions for analysis. For the 20% surface coverage models, the binding energy is relatively high, due to the AuAu bonds that are formed between the NP cores. This binding energy decreases with surface coverage and is only surpassed when the surface coverage reaches 80% (CH3) or more (NH2, OH). The increasing binding energy for high surface coverage NP model systems is due to the intra- and interparticle interactions of the surface ligands and the larger number of monomermonomer interactions. The AuAu interactions are minimized at about 50% surface coverage but are still appreciably lower at 80% and 100% surface coverage than at 20% surface coverage, except for the CH3 functionalized ligands due to corrugation. This indicates that fewer AuAu bonds are forming and passivation does occur as ligands are absorbed onto the Au surface. A slight increase in AuAu interactions above 80% coverage is observed in a breakout of the potential energy contributions and is due to corrugation of the Au surface. Corrugation was previously observed to be dependent upon surface coverage fraction.18 Binding Energy and Minimum Separation. One result of this effort has been the development of interaction potentials for use in mesoscale simulations. The first step in developing these interaction potentials is to compute the minimum packing distance and binding energy between gold NPs with different ligand surface densities and ligand terminal groups. To this end, we have investigated three terminal groups (CH3, OH, NH2) with
three different contact configurations and four surface coverage densities. Once the binding energies are computed, it is possible to compute derived values such as the Hildebrand solubility parameter. The NPNP binding energy and minimum separation distance were initially calculated using the two NPs in a vacuum model from the hysteresis analysis. From the simulation results shown in Figure 4, a number of observations can be made. First, for all ligand functional groups, a minimum binding energy occurs at 50% ligand surface coverage, when compared with the other surface coverages computed, namely, 20%, 80%, and 100%. This suggests that passivated gold NPs coated with alkanethiol ligands produce optimal interparticle separation when each particle is nearly 50% coated. Further simulations at similar coverage fractions are required to determine whether this is the minimum result. The binding distance does not appear to have a minimum value but rather increases monotonically from minimal surface coverage to 100% surface coverage. This is reasonable, since as surface coverage increases the thickness of the ligand coating also increases.18,31 The binding distance between the NPs is predictable from the polarity of the functional group as shown in the two-particle simulations. Specifically, the NP core separation is greatest for the most polar (OH) terminal group and tends to decrease with decreasing polarity (NH2 > CH3). Hildebrand Solubility Parameter. The Hildebrand solubility parameter is used to determine the miscibility behavior of nonpolar solvents. It can also be used to define the binding energy and core separation distance in a single number for use in mesoscale simulations. The Hildebrand solubility parameter is valid for the ligands and NP systems considered here, as the polarity of these ligands are a negligible portion of the system potential energy. The solubility parameter is defined as the square root of the cohesive energy density, eq 7. δ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðΔHvap RTÞ=Vm
ð7Þ
In eq 7, ΔHvap is the heat of vaporization, R is the ideal gas constant, T is temperature, and Vm is the molar volume. The molecular dynamics calculation used here, eq 8, is based on the cohesive energy density (CED) method described by Belmares et al.26 The simulations used to compute the solubility parameter are isothermalisobaric (NPT) ensembles. Unlike the previously discussed simulations, these simulations use periodic boundaries in order to compute the system pressure. The unit 7840
dx.doi.org/10.1021/la2005024 |Langmuir 2011, 27, 7836–7842
Langmuir
ARTICLE
Figure 5. Plot showing the computed solubility parameter results for each ligand functional group and versus surface coverage.
cell is therefore a well-defined volume for these calculations. In the previous simulations, an assumed BCC crystal structure was used to compute the packing density of the functionalized NPs. This assumption is important to the utilization of eq 8. Note that all of the following discussions of the solubility parameter use eq 8 to calculate values based on the molecular dynamics simulation results. v0 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 u n u E E h i i c C uB uBt ¼ 1 C C ð8Þ δ ¼ uB t@ Vc A NA n
∑
In eq 8, Ææ indicates a time average over the duration of the dynamics26 and Ei is the energy of individual NPs. Ec is the energy of a single unit cell, or system potential energy for the constant pressure and temperature periodic simulation. NA is Avogadro’s number, Vc is the unit cell volume, and n is the number of functionalized NPs. The Hildebrand solubility parameter is valid for nonpolar materials; this is the situation here where the alkanethiol ligands and the polymer solvents considered are primarily nonpolar. Therefore, this parameter should give a good estimate of the relative miscibility of the ligand-coated NPs. Equation 9 is the Hansen solubility parameter used elsewhere.26 Equation 9 is very similar to eq 8 except that k = 1,2,3 refers to the coulomb (polar), van der Waals (dispersion), and hydrogen bond contributions, respectively. For the polar terminal groups, namely, NH2 and OH, the Hansen solubility calculation, eq 9, is generally used. In this case, though the polar and hydrogen bond components make up less than 1% of the system energy, the Hildebrand solubility is therefore considered applicable. We have therefore chosen to ignore these component contributions for now, or in the case of the hydrogen bond the energy is already taken into account by the interatomic potential used. 0 n 1 k
k Ei Ec C B B C ð9Þ δk 2 ¼ Bt ¼ 1 C @ A Vc NA n
∑
In order to the compute the solubility of the alkanethiol ligandcoated gold NPs using the previous simulations, an assumption
must be made as to their packing arrangement. From Landman and Luedtke,10 we know that alkanethiol coated gold NPs may form a BCC nanocrystal array32 allowing for an estimate of the lattice binding energy. The cohesive energy of the superlattice per NP is 8.24ε, where ε is the binding energy between a pair of NPs.33 By assuming BCC packing, it is also possible to predict the corecore separation of the NPs and the binding energy from the NPT simulation results. After the cohesive energy density is computed, the solubility parameter is then computed as the square root of the cohesive energy density. In Figure 5, the Hildebrand solubility is plotted for each type of ligand functional group with respect to the surface coverage fractions. These plotted results come from the NPT simulations, but the results from the previous NVT simulations with the assumed BCC structure show similar trends and values. Solubility Results. In Figure 5, the Hildebrand solubility parameter results for each surface coverage and functional group are plotted. Additionally, a range of typically reported solubility values for polystyrene are also shown. A frequently reported range of solubility values for polystyrene is included in Figure 5, as it is a commonly used component of block copolymer matrices found in nanoparticle composites. Upon initial inspection, Figure 5 shows that the lowest computed solubility parameter for each functional group corresponds to the lowest computed binding energy in Figure 4b, of around 50%. This is an interesting result because it implies that it is possible to tune the NP solubility to a specific value by adjusting the surface coverage fraction. Having a tunable solubility will assist in developing a processing method for the manufacture of NP composites with a specific polymer matrix. Analysis of the observed minimum cohesive energy, or solubility parameter, has shown there to be competing interactions. The following discussion pertains primarily to CH3 and NH2 functional groups. The OH functional group results will be discussed following the current discussion. With minimally coated NPs (
