Molecular Simulation Studies on the Elongation of Gold Nanowires in

The George Washington University. , §. Present address: Division of Intramural Research, National Institute of Diabetes, Digestive, and Kidney Diseas...
0 downloads 0 Views 3MB Size
Download PDF
J. Phys. Chem. C 2010, 114, 10365–10372

10365

Molecular Simulation Studies on the Elongation of Gold Nanowires in Benzenedithiol Qing Pu,† Yongsheng Leng,*,‡ Xiongce Zhao,§,| and Peter T. Cummings†,| Department of Chemical and Biomolecular Engineering, Vanderbilt UniVersity, NashVille, Tennessee 37235, Department of Mechanical and Aerospace Engineering, The George Washington UniVersity, Washington, D.C. 20052, and Nanomaterials Theory Institute, Center for Nanophase Materials Sciences, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831 ReceiVed: February 24, 2010; ReVised Manuscript ReceiVed: May 2, 2010

The bonding geometry at the metal-molecule interface plays an important role in determining the conductance behavior of metal-molecule-metal junctions. This bonding geometry has to be determined a priori in quantum mechanical current-voltage (I-V) calculations. To identify the detailed metal-molecule bonding configurations, we perform classical molecular simulations by combining grand canonical Monte Carlo (GCMC) sampling with molecular dynamics (MD) to explore the dynamic elongations of gold nanowires in the presence of benzenedithiol (BDT) molecules. A specific multitime-scale double reversible reference system propagator algorithm (double-RESPA) has been designed for the metal-organic complex in MD simulations to improve the simulation efficiency. We investigate the variations of bonding sites and bonding angles of BDT molecules on a stretched Au nanowire at a constant chemical potential. The density of BDT and the number of bonded and nonbonded BDT molecules in the simulation box is monitored during the entire elongation process. Simulation results show that BDT molecules can form a denser monolayer on Au nanowires than at the Au (111) surface, owing to the many atomic steps on curved surfaces. Moreover, the chemical bonding of BDT on the Au nanowire significantly effect the elongation behavior of Au nanowires compared with those in vacuum. Our present results will be valuable to the understanding of the broken junction mechanism in molecular electronics conductance measurements. 1. Introduction Understanding molecule-metal contact effects on electron transport through metal-molecule-metal junctions is an important step toward building sophisticated molecular electronic devices.1 The transport properties of gold-thiolate-gold model systems have been studied experimentally in the past decade.2–4 In particular, the current-voltage (I-V) characteristics for benzenedithiolate (BDT) molecules sandwiched between two gold electrodes have been studied both experimentally3–6 and theoretically.7–11 One outstanding question concerns the discrepancy in I-V values between theory and experiments. Early first-principle calculations7 obtained more than 2 orders of magnitude in conductance larger than the experimental value,3 but comparable to other experimental results.6 It has been shown that the magnitude of the current is an extremely sensitive function of the contact geometry and chemistry in the metal-molecule interface. Temperature effect and local disorder of gold atoms all alter the value of the current.7 Clearly, the nanoscale contact geometry at the BDT-Au interface is one of the critical factors that have a dramatic influence on the conductance. For example, it has been found that conductance is dependent on whether the molecules are properly bonded to one or both gold electrodes at molecular ends.3 Binding sites of the thiol groups on the gold electrode and the molecular * Author to whom correspondence should be addressed. E-mail: [email protected]. † Vanderbilt University. ‡ The George Washington University. § Present address: Division of Intramural Research, National Institute of Diabetes, Digestive, and Kidney Diseases, National Institutes of Health, Bethesda, MD 20892. | Oak Ridge National Laboratory.

orientation with respect to the gold surface also influence the conductance.6,8 Several ab initio density functional theory (DFT) calculations12–14 showed different low energy configurations of the Au-BDT-Au system and stable binding sites. The S groups in BDT molecules can be at a bridge site, slightly displaced toward the fcc site or on a fcc sites. However, the difference in conductance due to different binding sites can be as large as 1 order of magnitude.13,14 In a mechanically controllable break junction (MCBJ) experiment, Tsutsui et al.15 discovered low conductance (∼0.01 G0) and high conductance (∼0.1 G0) regimes for BDT molecules bonded on Au wires. They attributed the two distinct conductance states to the Au-S bonding configurations changing from top-top to top-hollow (or the hollow-hollow) sites. This 10fold increase in conductance for Au-BDT-Au junction is consistent with Nara et al.13 theoretical calculations when the molecular bridging configuration changes from top-top to hollow-hollow sites. Using scanning tunneling microscope (STM) break-junction setup for single alkanedithiols, Li et al.16 also reported the occurrence of low and high conductance regimes. From a molecular perspective, it is valuable to perform classical molecular simulations to identify the detailed Au-BDT-Au configurations. This is because any I-V calibrations based on ab initio quantum mechanical calculations have to rely on realistic contact geometries at metal-molecule-metal interfaces. In this paper, we report our further molecular simulation studies on the dynamic elongations of gold nanowires in the presence of BDT molecules. We combine a grand canonical Monte Carlo (GCMC) simulation technique with molecular dynamics (MD) simulations to investigate the bonding structure and elongation behavior of gold nanowires in BDT solvent. A double reversible reference system propagator algorithm (double-

10.1021/jp101689u  2010 American Chemical Society Published on Web 05/25/2010

10366

J. Phys. Chem. C, Vol. 114, No. 23, 2010

RESPA) based on the original work of Tuckerman et al.17 and our previous work18 has been developed to increase the accuracy in multitime-scale MD simulations. The current work provides insights into the local bonding structures of BDTs on Au nanowires. Usually, DFT calculations are performed to search optimized geometries at 0 K, which are clearly subject to possible changes at ambient conditions in a complex environment. Our molecular simulation results can be used as an accurate input to quantum calculations that will help to resolve the discrepancies between theoretical and experimental predictions. Additionally, molecular modeling also helps to interpret and to give useful feedback to the experimental measurements, device design, and manufacturing. In section 2 we will first give a detailed description on the simulation methods, including GCMC and MD simulation techniques, force field parameters used in the simulations, and the double-RESPA. Section 3 will present our simulation results and discussion, followed by conclusions in section 4. 2. Simulation Methods 2.1. Grand Canonical Monte Carlo (GCMC) Simulation. To simulate the adsorption processes of BDT molecules on gold nanowires, the GCMC sampling technique is used. The BDT molecules are modeled as rigid molecules from an optimized structure using universal force field (UFF).19 GCMC moves involve the insertion and deletion of BDT molecules for a given chemical potential (µ), as well as the translation and rotation of rigid BDT molecules according to a set of given probabilities. Simulations are carried out at 298 K. The simulation box lengths in the x- and y-dimensions are 3.5 and 3.0 nm, respectively, with a varying box length in the z-direction, or the [001] elongation direction. Periodic boundary conditions are applied in all three directions. Detailed descriptions on GCMC µVT sampling and subsequent canonical NVT sampling is given in our previous paper.20 2.2. Combining GCMC and Molecular Dynamics (MD) Simulations. For the elongation of an Au nanowire along the [001] direction in BDT solvent, we propose a method of combining GCMC and MD techniques alternatively to keep the BDT density at a bulk value during the elongation process. This is because in the original experiment by Reed et al.,3 elongation of the gold junction was performed through mechanical controllable break junction experiment in BDT solutions.21 To simulate this process in the fabrication environment, initially, the GCMC simulation is performed to obtain an equilibrium adsorption structure of BDT on a gold nanowire (see Figure 2a). This equilibrium structure is then taken as the initial structure for the MD simulations of the elongation of Au nanowire in BDT, typically with a small pulling distance of 5 Å along the [001] direction. The elongation of the nanowire is realized by moving the top rigid layers along the z-direction in increments of 0.1 Å, followed by 1000 time step relaxations, with a time step of 2 fs. This corresponds to an elongation rate of 5.0 m/s. The temperature is controlled at 298 K. During the elongation process, we assume that the binding sites of BDT molecules on an Au nanowire are unchanged. After the 5 Å elongation MD simulation is finished, MD is switched back to GCMC for a new equilibrium adsorption run. The resulting adsorption configuration is then used for a new MD elongation run. These two procedures proceed alternatively until the breakage of the gold nanowire in the BDT solvent. Additional BDT molecules in GCMC runs will be inserted into the simulation box to keep the BDT density equal to the bulk density. According to the experimental measurements performed by Reed et al.,3 after the breakage of the gold nanowire, the solvent was allowed to

Pu et al. evaporate under ambient conditions. We also attempt to mimic this process by removing all nonbonded BDT molecules before bringing the tips together. The formation mechanism of this Au-BDT-Au junction during the bringing in process is still under investigation. 2.3. Force Field Parameters in Molecular Simulations. In the MD simulation part of elongation and breakage of a gold nanowire in BDT solvent, BDT molecules are modeled as flexible molecules. Following our previous study,22 the total potential energy of metal-organic complex is expressed as a sum of several energy terms including chemical bonded interactions within BDT molecules, and between Au and BDT, as well as nonbonded interactions among BDT and Au interactions. The bonded interactions within BDT consist of bond stretching, bond angle bending, dihedral angle torsion, and inversion terms. The nonbonded interactions consist of van der Waals and electrostatic terms. A 12-6 Lennard-Jones potential applies to atoms in the same BDT molecule but separated by more than two neighboring atoms (the 1-2 and 1-3 interaction exclusions) or for atoms in different BDT molecules, as well as for all of the nonbonded Au atoms interacting with BDT molecules. The Mulliken partial charges from ab initio quantum mechanical calculations for BDT molecules are also considered in electrostatic energy calculations.22 The chemical bonded interactions between the S atom in BDT and its bonded Au atoms are described by a Morse potential. On the basis of our previous work,22 the Morse potential takes the form Au-S UMorse ) E0 exp[-R(r - r0)]{exp[-R(r - r0)] - 2}

(1) where r is the distance between the S atom and its bonded Au atoms. The force parameters E0, r0, and R are fitted from DFT energy calculations.22 The Morse potential parameters have been further refined recently on the basis of the DFT calculations for different BDT-nAu complexes.23 The on-top and on-bridge bonding sites are two cases considered in our present work. In particular, the Morse potential parameters for BDT-3Au12 and BDT-3Au21 configurations23 are used for the on-top and onbridge bonded interactions, respectively. These force parameters are listed in Table 1. The BDT-Au angle bending and rotation potentials are not considered in our present study, since they are largely dominated by the intermolecular interactions.22 In dealing with the occurrence of chemical bonding between S atoms in BDT molecules and surface gold atoms, we assume that the interaction depends on the S-Au distance. If neither of the S atoms in a BDT molecule is within the Morse interaction distance (r < 2.5r0) described in eq 1 with Au, the S-Au interaction is calculated by the pairwise Lennard-Jones potential. As soon as one of the S atoms in a BDT molecule is in the range of 2.5r0 from any of the surface Au atoms, the Morse potential is in effect. Simultaneously, the Lennard-Jones interaction is turned off, and the H atom originally bonded to the S atom will be removed from the whole simulation system. In previous molecular simulation studies,18,22 we demonstrated that the S-H bond in BDT is quite flexible at room temperature (the rotation barrier from perpendicular to parallel configuration to benzene ring is comparable to kT). This S-H bond has almost minimal effect on the rigid planar structure of BDT molecules. The Au-Au atomic interactions are described by the secondmomentum approximation of tight-binding (TB-SMA) potential.24 We have recently tested several atomic potentials for

Elongation of Gold Nanowires in Benzenedithiol

J. Phys. Chem. C, Vol. 114, No. 23, 2010 10367 adsorbed by the external bath during the equilibration. Mathematically, any integrators derived from the Trotter factorization of the Liouville operator iL are reversible.17 On the basis of our previous work18 for disparate mass problems, we derived a double-RESPA algorithm for the Au-BDT interaction system, which involves a further decomposition of the propagator for H atoms in BDT and separations of positions and momenta between BDT particles and Au atoms. The detailed numerical procedure of the double-RESPA is presented in the Appendix.

Figure 1. Variations of different energy terms in double-RESPA for the equilibration of BDT molecules on a 256-atom Au (001) nanowire during the first 10 000 time step MD run with ∆t ) 2 fs.

TABLE 1: Force Parameters in the Au-S Morse Potential23 E0 (kcal/mol)

r0 (Å)

R (1/Å)

28.7 49.39

2.44 2.29

1.67 1.72

bridge site top site

transition metals and concluded that TB-SMA is the most suitable one for the description of energetics for small gold clusters.25 2.4. Double Reversible Reference System Propagator Algorithm (Double-RESPA). Due to the high stiffness of the benzene ring and strong valence bonding in BDT molecules, in MD simulation part a much smaller time step is needed to accurately integrate the fast dynamics of BDT contributed from strong intramolecular forces, while at the same time a comparably larger time step is enough to integrate the slow dynamics contributed from slowly varying intermolecular forces. Moreover, within the intramolecular forces, the H atoms are 1 order of magnitude smaller in mass than other particles, resulting in even higher vibration frequencies than others. Therefore, a multiple-time-scale method involving BDT and heavy gold atoms with different time step integrations is appropriate to integrate the equations of motion of the whole system. In our present study, we have further reformulated the double reversible reference system propagator algorithm (RESPA) for Au-BDT complex. It is based on the original work of Tuckerman et al.17 and our previous work,18 in the extended Nose´-Hoover dynamics26,27 for the NVT ensemble. The conserved quantity (pseudo-Hamiltonian) of the system is given by

H'(x,y,px,py,η,pη) ) V(x,y) +

∑

px2 + 2m

∑

py2 pη2 + + 3NkTη ) V + 2m 2Q Vkin + VNose (2)

where x and y are the positions of BDT particles and Au atoms, respectively, and px and py are their related momenta. Variables η and pη are the extended thermostat and conjugated momentum. Quantities N, k, and T are the total number of MD particles, the Boltzmann constant, and temperature. In eq 2, V is the total potential energy of the molecular system, and Vkin is the kinetic energy. The last term, VNose, can be recognized as the heat

3. Results and Discussion 3.1. Double-RESPA Validation. We first check the validity of the double-RESPA scheme. Figure 1 shows the variations of different energy terms in eq 2 during the first 20 ps MD relaxation after the 0.1 Å pulling of Au nanowire. The time step in the outer-loop is ∆t ) 2 fs. In implementing the doubleRESPA scheme, the numbers of inner and inner-inner loops are n ) 5 and m ) 4; therefore, the smallest time step for H atoms is 1/20th of the outer-loop time step, i.e., 0.1 fs. The mass of the thermostat Q is chosen according to the relation Q ) 3NkTτ, where τ is a characteristic time in the system.17 We find that MD equilibration is not sensitive to the value of Q. Therefore, we arbitrarily choose τ ) 500 fs in this study. From Figure 1 we see that the system approaches equilibrium in about 1000 time steps (i.e., 2 ps). The pseudo-Hamiltonian, H′, remains constant except for a small fluctuation in the first few time steps. Therefore, the modified double-RESPA scheme works well for our Au-BDT system. 3.2. BDT Bonding Configurations in Elongation Process. To understand how the bonding configuration of BDT molecules on a gold nanowire changes during the elongation process, we investigate the variations of bonding sites of BDT on the gold nanowire, the BDT-Au bonding angle and other geometric parameters versus elongation distance, while keeping the bulk density of BDT constant through GCMC insertion and deletion of BDT molecules. Initially, a GCMC run with ten million (10 M) moves is performed, beginning from an artificial configuration in which a bare 256-atom Au (001) nanowire is surrounded by a few BDT molecules. We monitor the change of the number of BDT molecules during the GCMC insertions and deletions for a given chemical potential, until it reaches a stable value. After 10 M GCMC moves are finished, another canonical MC run with 10 M moves is carried out for further equilibration of the system. These two procedures are consequently applied upon every MD elongation with 5 Å increment is finished. 3.2.1. BDT Bonding Sites. Figure 2 shows the side and top views of the initial equilibrium adsorption structure of BDT molecules on a 256-atom Au (001) nanowire before elongation. The total number of BDT molecules in the simulation box at initial equilibrium is 142, of which 60 BDT molecules are found to be bonded on the gold nanowire. We find that only one ontop bonding site is formed, the rest are all on-bridge bonding. Table 2 lists the calculated BDT density in the simulation box after each 5 Å increment elongation. Obviously, the BDT density remains the same bulk value during the entire elongation process of the gold nanowire, reflecting the real situation in experimental measurements.21 We find that the total number of BDT molecules (including both bonded BDT molecules on the Au nanowire and nonbonded BDT molecules in the bulk) increases by approximately 25 after each 5 Å incremental elongation and subsequent 10 M GCMC runs. Of these 25 added BDT molecules, the portion of newly bonded BDT molecules exhibits a deceasing trend. That is, the number of newly bonded BDTs decreases from 15 at 5 Å elongated state to 5 at 20 Å elongated

10368

J. Phys. Chem. C, Vol. 114, No. 23, 2010

Pu et al. TABLE 3: Distribution of BDT Molecules Bonded on the 256-Atom Au Nanowire after Each 5 Å Incremental Elongation along the [001] z-Direction

Figure 2. Side and top views of the initial equilibrium adsorption configuration of 142 BDT molecules on a 256-atom Au (001) nanowire prior to pulling. The dimension of the simulation box is 3.5 × 3 × 3.3 nm3.

TABLE 2: BDT Densities in Simulation Box during the Elongation of the 256-Atom Au Nanowire along the [001] z-Direction elongation (Å)

BDT density (g/cm3)

5 10 15 20

1.06 1.10 1.10 1.10

state, as indicated in Table 3. In other words, most newly added BDT molecules in the simulation box are nonbonded bulk molecules. This indicates that the gold nanowire under elongation, particularly in the much elongated stages, will not generate more surface atoms. The stretched neck observed in simulations only accommodates a few more bonded BDT molecules. Table 3 also shows that there are more on-bridge bonding sites than on-top bonding sites, suggesting that on-bridge bonding of BDT on gold nanowires is more favorable. This is particularly true for the Au nanowire before elongation. Similarly to our previous studies,22,23,28 we calculate the distance distribution functions of S-Au and S-S on the Au nanowire. The two distribution curves are shown in Figure 3. The first S-Au peak (Figure 3a) seems always at a universal distance of 2.5 Å,22 regardless of whether it is on-top or onbridge bonding. The second and third S-Au peaks at larger distances correspond to the nonbonded gold surface atoms. The first S-S peak shown in Figure 3b still occurs at about 4.0 Å. This S-S distance is much smaller than the 5.0 Å distance commonly observed for thiol molecules on a Au (111) surface.29–31 Most likely, this is attributed to the many atomic steps on the curved surface of the gold nanowire that favors a denser BDT

elongation (Å)

on-top bonding

on-bridge bonding

total bonded BDT

0 5 10 15 20

1 30 30 34 36

59 45 55 56 59

60 75 85 90 95

monolayer bonded on the surface (see Figure 2). Interestingly, Figure 3b further shows that as the gold nanowire is elongated substantially and becomes thinner, the S-S distribution curve shifts toward an even smaller distance, indicating that more bonded BDT molecules aggregate on the thinner Au nanowire. For S-Au distribution (Figure 3a), the first peak of S-Au distance always stays at 2.5 Å. To see how the BDT adsorption structure evolves with the elongation of the gold nanowire, Figure 4 shows the three snapshots of BDT on the 256-atom gold nanowire at different elongation stages at room temperature. 3.2.2. BDT-Au Bonding Angle. The BDT-Au bonding angle determines the electron density between S and bonded Au atoms and is closely related to the electron transport properties of the metal-molecule-metal junction if the molecule is connected to two gold electrodes.8,13,14,32–35 Usually, the orientation of thiolate molecules on the Au (111) surface of the gold electrode is described by the tilt angle. For very small gold nanowires, direct determination of the surface normal direction is difficult. Here we turn to the calculation of the C-S-Au angle. In general, the tilt angle is related to, but not exactly equal to 180° - ∠C-S-Au. To understand how the adsorption configuration changes during the entire elongation process, we investigate the probability distribution of the angle ∠C-S-Au in two situations. One is at the initial equilibrium configuration prior to pulling (Figure 2a) and the other at intermediate elongation in which a thinner neck of gold nanowire is formed after 15 Å pulling (Figure 4c). Statistical averages are obtained over a 100 ps MD run after the system reaches equilibration. The results of the probability distributions of angles have been separated into on-top and on-bridge bonding cases. Figure 5 shows that there is a wide spectrum of angle distributions, regardless of on-top or on-bridge bonding. Taking the probability density at 0.01 as a reference line, we see that a large portion of C-S-Au angles before pulling are located between 110-140° for on-top bonding and 100-150° for on-bridge bonding. However, at the elongated stage (15 Å elongation), these values are decreased by 20-30°. This

Figure 3. S-Au (a) and S-S (b) distance distribution functions at the initial and intermediate elongation stages.

Elongation of Gold Nanowires in Benzenedithiol

J. Phys. Chem. C, Vol. 114, No. 23, 2010 10369

Figure 4. Snapshots of BDT on a 256-atom Au nanowire elongated in bulk BDT solvent at room temperature: (a) ∆z ) 5 Å; (b) ∆z ) 10 Å; (c) ∆z ) 15 Å.

Figure 5. Angle C-S-Au density distributions at (a) on-top bonding and (b) on-bridge bonding of BDT on Au nanowire for initial configuration prior to pulling and for intermediate configuration after 15 Å elongation.

Figure 6. Snapshots of BDT on a 256-atom Au nanowire broken in bulk BDT solvent and subsequent bringing to form the Au-BDT-Au junction at room temperature: (a) breaking at ∆z ) 23 Å; (b) reconnecting at ∆z ) -5 Å; (c) reconnecting at ∆z ) -10 Å.

decreasing in the C-S-Au bonding angle is most likely due to the stronger BDT-BDT intermolecular interactions at very high-curvature surfaces that have many defects on the surfaces.

3.3. Effect of BDT on the Ductile Elongation of Gold Nanowires. The Break-Junction Structure and Forces. In our previous work on the elongation of gold nanowire in nonbonded propane solvent, we found that the solvent

10370

J. Phys. Chem. C, Vol. 114, No. 23, 2010

molecules have almost no effect on the ductile elongation of gold nanowire. When the gold nanowires are immersed in BDT, the ductile elongation behavior has been found significantly depressed in the present simulation study. In the vacuum at an elongation rate of 5.0 m/s, we obtained a 15-atom-long monatomic chain structure at 298 K.36 The total ductile elongation of the gold nanowire is 45.4 Å. In BDT solvent, we find that the ductile elongation is decreased to 23 Å in this specific case. Clearly, the surrounding BDT molecules, particularly the bonded BDT molecules, have a significant impact on the breakage of Au nanowires. Occasionally, we found that a small gold cluster with a few gold atoms was pulled out with bonded BDT molecules during the breakage, indicating that the S-Au bond is stronger than the Au-Au bond in this specific situation. This observation supports recent density functional based molecular dynamics simulation results,37,38 as well as experimental predictions.39 Figure 6 shows the three snapshots of BDT molecules on a 256-atom Au nanowire at break, and at the subsequent bringing to contact to form the Au-BDT-Au junction at room temperature. According to the experimental measurements performed by Reed et al.,3 after the breakage of the gold nanowire, the solvent was allowed to evaporate under ambient conditions. We mimic this process by removing all nonbonded BDT molecules before bringing the BDT monolayer covered tips together (see Figure 6b,c). The formation mechanism of Au-BDT-Au junction at this stage is still under investigation. We also calculate the breaking forces of the gold nanowire in both vacuum and BDT solvent at room temperature. Force calculations are based on the averaged values in the final MD equilibrium run before break. Due to the thermal fluctuations, the breaking forces of the gold nanowire, either in vacuum or in BDT, are around 0.5 nN, which is much smaller than what we obtained at low temperature (∼1.5 nN).25 4. Conclusions Combining grand canonical Monte Carlo (GCMC) simulation with multiple-time-scale MD technique (double-RESPA), we explore the elongation behavior of gold nanowires in BDT solvent and the related adsorption structure. BDT molecules form chemical bonds with Au atoms, either on the on-top or on-bridge site. Simulations of Au nanowire elongations are performed at the constant chemical potential (or at the BDT bulk density), reflecting the real situation of stretching gold junctions in experiment. The large curvature surface and many surface steps on the stretched Au nanowire promote the formation of a denser BDT monolayer on the surface, indicated by the shifts of the C-S-Au angle and S-S distribution curves during the elongation of the gold nanowire. It seems that the effect of bonded BDT molecules on gold nanowires results in an earlier breakage of gold nanowires than in vacuum. Acknowledgment. This work was supported by the U.S. Department of Energy (DOE) Office of Science the Computational Nanoscience project, and also by the National Energy Research Scientific Computing Center (NERSC). Appendix: Double-RESPA The double reversible reference system propagator algorithm (double-RESPA) is based on the original work of Tuckerman et al.17 and has been described in detail in our previous work.18

Pu et al. For the BDT-Au system, we need to rewrite the Liouville operatorcorrespondingtotheextendedNose´-Hooverdynamics26,27 (i.e., eq 2) as the following:17

[

iL ) x˙

]

∂ ∂ ∂ - η˙ px + + Fx(x,y) ∂x ∂px ∂px ∂ ∂ ∂ - η˙ py + y˙ + Fy(x,y) ∂y ∂py ∂py ∂ ∂ η˙ ) iLr + + Fη(p) ∂η ∂pη ∂ [Fl(x,y) - η˙ px] + ∂px ∂ ∂ y˙ + + [Fy(x,y) - η˙ py] ∂y ∂py ∂ ∂ η˙ ) iLr + iL1 + Fη(p) ∂η ∂pη

[

]

(A1)

where

Fη(p) ) p˙η )

2

∑ pm - 3NkT

and coordinates x and y refer to the positions of atoms in BDT and Au atoms, respectively, and px and py are their related momenta. Variables η and pη are the extended thermostat and conjugated momentum. iLr is the reference system operator and iL1 is the remaining operator. iLr is defined as

iLr ) x˙

∂ ∂ + Fs(x) ∂x ∂px

(A2)

where Fs(x) is the short-range intramolecular interactions within BDT molecules that are only dependent on the local coordinates, x, of BDT molecules. In eq A1, Fl(x,y) ) Fx(x,y) - Fs(x), is the long-range intermolecular interactions between BDT molecules, as well as between BDT and Au atoms. Suppose that Γ(0) ) Γ(x(0), y(0), px(0), py(0), η(0), pη(0)) is the initial state of the molecular system. The state of the system at time t is given by

Γ(t) ) eiLtΓ(0)

(A3)

For a small time step ∆t, the propagator eiL∆t can be factored using the Trotter theorem as17

eiL∆t ) eiL1∆t/2eiLr∆teiL1∆t/2

(A4)

The middle propagator in eq A4 involves the calculation of short-range forces that have fast varying characteristics with large magnitudes. This will need a smaller time step δt ) ∆t/n to integrate the equations of motion in the reference system for n times. Because of the large difference in mass between H atoms and C, S atoms, to avoid the drift in total energy, we follow the same procedure18 to further decompose the propagator for light particles and heavy particles in BDT within the reference system:

iLr ) iLlr + iLhr

(A5)

Elongation of Gold Nanowires in Benzenedithiol

J. Phys. Chem. C, Vol. 114, No. 23, 2010 10371

Y ) e∆t/2[y˙(∂/∂py)+Fy(x,y)(∂/∂py)-η˙ py(∂/∂py)] ) e-∆t/4η˙ py∂/∂pye∆t/4Fy(x,y)∂/∂pye∆t/2y˙∂/∂ye∆t/4Fy(x,y)∂/∂pye-∆t/4η˙ py∂/∂py (A13)

where

∂ ∂ iLlr ) ˙l + Fls(l,h) ∂l ∂pl

∂ iLhr ) h˙ + ∂h

and

Fhs(l,h)

∂ ∂ph

Here, l and h label light and heavy particle coordinates. The middle propagator can then be written as

eiLr∆t ) [eiLlr(δt/2)eiLhδteiLlr(δt/2)]n

(∆t ) nδt)

(A6) In the above equation, a second decomposition is employed in the propagator eiLlr(δt/2), yielding

eiLlr(δt/2)eiLhrδteiLlr(δt/2) ) [eiLlr(dt)]m/2eiLhrδt[eiLlr(dt)]m/2 (δt ) mdt) (A7) where ˙

eiLlr(dt) ) e(dt/2)Fls(l,h)∂/∂pledtl∂/∂le(dt/2)Fls(l,h)∂/∂pl ˙

eiLhrδt ) e(δt/2)Fhs∂/∂pheδth∂/∂he(δt/2)Fhs∂/∂ph

(A8)

(A9)

The total gain of this double-RESPA scheme is mn. The detailed propagator algorithms within the reference system, eiLr∆t (the middle propagator in eq A4), are given in ref 18. The key issue in developing double-RESPA based on the propagator (A4) is to appropriately decompose the side propagator eiL1∆t/2. We rewrite the remaining Liouville operator from (A1) as

iL1 ) [Fl(x,y) - η˙ px]

∂ ∂ ∂ + y˙ + + [Fy(x,y) - η˙ py] ∂px ∂y ∂py ∂ ∂ η˙ (A10) + Fη(p) ∂η ∂pη

There are many factorization algorithms for eiL1∆t/2 that satisfy the Trotter theorem.17 Our choice of factorization is to properly update the positions and forces of BDT and Au atoms evenly, and at the same time to keep the total pseudo-Hamiltonian, H′, at a constant energy surface. The final decomposition of the propagator (A4) is given as

eiL∆t ) e∆t/2Fη(p)∂/∂pηXYXe∆t/2η˙ ∂/∂ηeiLr∆te∆t/2η˙ ∂/∂ηXYXe∆t/2Fη(p)∂/∂pη (A11) where

X ) e∆t/4[Fl(x,y)(∂/∂px)-η˙ px(x,y)(∂/∂px)] ) e∆t/8Fl(x,y)∂/∂pxe-∆t/4η˙ px∂/∂pxe∆t/8Fl(x,y)∂/∂px and

Given the reference propagator algorithm in eqs A6-A9, and the sequence of all the propagators shown in eqs A11-A13, the detailed propagation of the phase space from initial Γ(0) ) {xl(0), xh(0), y(0), pxl(0), pxh(0), py(0), η(0), pη(0)} to Γ(∆t) ) {xl(∆t), xh(∆t), y(∆t), pxl(∆t), pxh(∆t), py(∆t), η(∆t), pη(∆t)} at time ∆t, will follow the routine procedures, as described in our previous work.18

(A12)

References and Notes (1) Nitzan, A.; Ratner, M. A. Science 2003, 300, 1384–1389. (2) Cui, X. D.; Primak, A.; Zarate, X.; Tomfohr, J.; Sankey, O. F.; Moore, A. L.; Moore, T. A.; Gust, D.; Harris, G.; Lindsay, S. M. Science 2001, 294, 571–574. (3) Reed, M. A.; Zhou, C.; Muller, C. J.; Burgin, T. P.; Tour, J. M. Science (Washington, D.C.) 1997, 278, 252–254. (4) Xu, B. Q.; Tao, N. J. J. Science 2003, 301, 1221–1223. (5) Tsutsui, M.; Shoji, K.; Morimoto, K.; Taniguchi, M.; Kawai, T. Appl. Phys. Lett. 2008, 92, 223110. (6) Xiao, X. Y.; Xu, B. Q.; Tao, N. J. Nano Lett. 2004, 4, 267–271. (7) Di Ventra, M.; Pantelides, S. T.; Lang, N. D. Phys. ReV. Lett. 2000, 84, 979–982. (8) Kornilovitch, P. E.; Bratkovsky, A. M. Phys. ReV. B 2001, 64, 195413. (9) Krstic, P. S.; Dean, D. J.; Zhang, X. G.; Keffer, D.; Leng, Y. S.; Cummings, P. T.; Wells, J. C. Comput. Mater. Sci. 2003, 28, 321–341. (10) Seminario, J. M.; Zacarias, A. G.; Tour, J. M. J. Phys. Chem. A 1999, 103, 7883–7887. (11) Stokbro, K.; Taylor, J.; Brandbyge, M.; Mozos, J. L.; Ordejon, P. Comput. Mater. Sci. 2003, 27, 151–160. (12) Kondo, H.; Nara, J.; Kin, H.; Ohno, T. Jpn. J. Appl. Phys. 2008, 47, 4792–4798. (13) Nara, J.; Geng, W. T.; Kino, H.; Kobayashi, N.; Ohno, T. J. Chem. Phys. 2004, 121, 6485–6492. (14) Pontes, R. B.; Novaes, F. D.; Fazzio, A.; da Silva, A. J. R. J. Am. Chem. Soc. 2006, 128, 8996–8997. (15) Tsutsui, M.; Teramae, Y.; Kurokawa, S.; Sakai, A. Appl. Phys. Lett. 2006, 89, 163111. (16) Li, X. L.; He, J.; Hihath, J.; Xu, B. Q.; Lindsay, S. M.; Tao, N. J. J. Am. Chem. Soc. 2006, 128, 2135–2141. (17) Tuckerman, M.; Berne, B. J.; Martyna, G. J. J. Chem. Phys. 1992, 97, 1990–2001. (18) Leng, Y. S.; Keffer, D. J.; Cummings, P. T. J. Phys. Chem. B 2003, 107, 11940–11950. (19) Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A.; Skiff, W. M. J. Am. Chem. Soc. 1992, 114, 10024–10035. (20) Zhao, X. C.; Leng, Y. S.; Cummings, P. T. Langmuir 2006, 22, 4116–4124. (21) Reed, M. A.; Zhou, C.; Muller, C. J.; Burgin, T. P.; Tour, J. M. Science 1997, 278, 252–254. (22) Leng, Y. S.; Krstic, P. S.; Wells, J. C.; Cummings, P. T.; Dean, D. J. J. Chem. Phys. 2005, 122, 244721. (23) Leng, Y. S.; Dyer, P. J.; Krstic, P. S.; Harrison, R. J.; Cummings, P. T. Mol. Phys. 2007, 105, 293–300. (24) Cleri, F.; Rosato, V. Phys. ReV. B 1993, 48, 22–33. (25) Pu, Q.; Leng, Y. S.; Tsetseris, L.; Park, H. S.; Pantelides, S. T.; Cummings, P. T. J. Chem. Phys. 2007, 126, 144707. (26) Hoover, W. G. Phys. ReV. A 1985, 31, 1695–1697. (27) Nose, S. J. Chem. Phys. 1984, 81, 511–519. (28) Pu, Q.; Leng, Y. S.; Zhao, X. C.; Cummings, P. T. Nanotechnology 2007, 18. (29) Nuzzo, R. G.; Allara, D. L. J. Am. Chem. Soc. 1983, 105, 4481– 4483. (30) Ulman, A. An Introduction to Ultra-thin Organic Films: From Langmuir-Blodgett to Self-Assembly; Academic Press: San Diego, CA, 1991. (31) Hautman, J.; Klein, M. L. J. Chem. Phys. 1989, 91, 4994–5001. (32) Haiss, W.; Wang, C. S.; Grace, I.; Batsanov, A. S.; Schiffrin, D. J.; Higgins, S. J.; Bryce, M. R.; Lambert, C. J.; Nichols, R. J. Nat. Mater. 2006, 5, 995–1002. (33) Nara, J.; Higai, S.; Morikawa, Y.; Ohno, T. J. Chem. Phys. 2004, 120, 6705–6711. (34) Letardi, S.; Cleri, F. J. Chem. Phys. 2004, 120, 10062–10068.

10372

J. Phys. Chem. C, Vol. 114, No. 23, 2010

(35) Kruger, D.; Fuchs, H.; Rousseau, R.; Marx, D.; Parrinello, M. J. Chem. Phys. 2001, 115, 4776–4786. (36) Pu, Q.; Leng, Y. S.; Cummings, P. T. J. Am. Chem. Soc. 2008, 130, 17907–17912. (37) Kruger, D.; Fuchs, H.; Rousseau, R.; Marx, D.; Parrinello, M. Phys. ReV. Lett. 2002, 89, 186402.

Pu et al. (38) Paulsson, M.; Krag, C.; Frederiksen, T.; Brandbyge, M. Nano Lett. 2009, 9, 117–121. (39) Huang, Z. F.; Chen, F.; Bennett, P. A.; Tao, N. J. J. Am. Chem. Soc. 2007, 129, 13225–13231.

JP101689U