3990
Langmuir 1997, 13, 3990-4002
Computer Simulation of the Adsorption of Alkanethiols on Au(111) from the Gas Phase. 1. Methanethiol H. Morgner Institut fu¨ r Experimentalphysik, Universita¨ t Witten-Herdecke, Stockumer Strasse 10, D-58448 Witten, Federal Republic of Germany Received October 28, 1996. In Final Form: February 12, 1997X Adsorption of methanethiol molecules from gas phase onto a Au(111) substrate is investigated by means of a molecular dynamics (MD) computer simulation. The problem that some of the elementary processes during adsorption take place on a time scale far too long to be handled by standard MD techniques is solved by making use of so-called “constrained or targeted dynamics”. It is found that the lifetime against desorption out of the chemisorbed precursor state depends strongly on the coverage. It varies from 38 to 0.6 µs between low and high coverage. The sticking coefficient is found to drop from >2 × 10-3 to below 1 × 10-5. By comparison with experimental data it is possible to determine a rate constant k ∼ 65 s-1 for the transition from the precursor state into the reactively adsorbed thiolate state.
1. Introduction Self-assembled monolayers (SAM’s) of organic molecules are a way to modify properties of solid surfaces.1 Alkanethiols HS(CH2)n-1CH3 adsorbed on metal surfaces represent a widely investigated prototype of such systems. The structure of alkanethiol films has been investigated by a variety of experimental techniques like electron diffraction,2 X-ray diffraction,3 and helium atom diffraction4 as well as optical ellipsometry, infrared spectroscopy, and electrochemistry.5 In most studies the substrate is a gold monocrystal surface, but other metals like silver6 and platinum7 have been used as well. The above quoted literature concentrates on the structural properties of completed monolayers. So far only little effort has been devoted to the kinetics of SAM formation. Buck et al.8 have followed the adsorption of an alkanethiol film from solution by means of optical second harmonic generation. This method is distinguished in that it allows in-situ observation. A comparison with ex-situ X-ray photoelectron spectroscopy (XPS) data gave satisfying agreement. As a result the authors state that the monolayer formation follows Langmuir adsorption kinetics. In the case of a system as complex as a monolayer formed of long chain alkanethiols it is doubtful that the measured kinetic behavior allows an unambiguous conclusion as for the adsorption mechanism. The fact that the observed order within a completed SAM is due to the interaction between neighboring adsorbates necessarily renders one of the basic assumptions of Langmuir kinetics invalid, namely, that the ability of a molecule to adsorb at a given site is independent of the occupation of neighboring sites.9 Indeed, another study points to a distinct non-Langmuir behavior. Bain et al.10 observed X
Abstract published in Advance ACS Abstracts, June 1, 1997.
(1) Ulman, A. An Introduction to Ultrathin Organic Films; Academic Press: Boston, MA, 1991. (2) Strong, L.; Whitesides, G. M. Langmuir 1988, 4, 546-558. (3) Fenter, P.; Eisenberger, P.; Liang, K. S. Phys. Rev. Lett. 1993, 70, 2447-2450. (4) Camillone, N., III; Chidsey, Ch. E. D.; Liu, G.; Putvinski, T. M.; Scoles, G. J. Chem. Phys. 1991, 94, 8493-8502. Camillone, N., III; Chidsey, Ch. E. D.; Eisenberger, P.; Fenter, P.; Li, J.; Liang, K.S.; Liu, G.; Scoles, G. J. Chem. Phys. 1993, 99, 744-747. (5) Porter, M. D., Bright, Th. B.; Allara, D. L.; Chidsey, Ch. E. D. J. Am. Chem. Soc. 1987, 109, 3559-3568. (6) Dubois, L. H.; Zegarski, B. R.; Nuzzo, R. G. J. Chem. Phys. 1993, 98, 678-689. (7) Bru¨ckner, M.; Heinz, B.; Morgner, H. Surf. Sci. 1994, 319, 37080. (8) Buck, M.; Eisert, F.; Fischer, J.; Grunze, M.; Tra¨ger, F. Appl. Phys. 1991, A53, 552-6.
S0743-7463(96)01043-8 CCC: $14.00
a fast accumulation of alkanethiols out of an ethanol solution onto gold up to 80-90% of maximum thickness, followed by a period of very slow adsorption which lasts several hours until the layer reaches its final thickness. Processes with very long characteristic times have been observed even after removal of alkanethiol films from solution. Poirier et al.11 observed by means of ultrahigh vacuum scanning tunneling microscopy (UHV-STM) that short chain homologues (n ) 4 and n ) 6) exhibit a 2D liquid phase which partly transforms into ordered domains having a unit cell of p×x3 (8 e p e 10) over periods of several hours. There is some other experimental information that mustsin principlescontain a clue as to adsorption mechanisms. The attempt to build up monolayers of alkanethiols on Au(111) by adsorption in UHV from the gas phase has been successful for methanethiol and alkanethiols with longer hydrocarbon chains up to decanethiol.6 In this context the “successful” film formation implies not only the deposition of molecules on the surface but the formation of a film with sufficient degree of lateral order to produce a recognizable pattern in low-energy electron diffraction (LEED). This criterion could not be fulfilled for alkanethiols HS (CH2)n-1CH3 with n > 10. Within the homologous series with n ) 1 ... 10 the authors report two interesting observations. The highest exposure to make an ordered film is needed for methanethiol (>104 langmuir). The necessary exposure drops with increasing chain length to reach a minimum of ∼10 langmuir with hexanethiol. Beyond n ) 6 again higher gas doses are required, the exposure needed for decanethiol being several thousand langmuirs. The observed LEED patterns point out the familiar x3×x3 structure for methanethiol and ethanethiol, but for the longer alkane thiols the authors report mx3×x3 unit cells, m ranging between 4 and 6. For the 5x3×x3 pattern the authors propose a unit cell structure which has only 80% of the x3×x3 monolayer coverage.6 In adsorption from liquid solution the shorter alkanethiols HS(CH2)n-1CH3 (n e 10) form the x3×x3 structure.12 Consequently, the mx3×x3 structure observed in gas phase adsorption does not represent the thermodynamic equilibrium but must be (9) Atkins, P. W. Physical Chemistry; Oxford University Press: Oxford, 1994. (10) Bain, C. D.; Troughton, E. B.; Tao, Y.-T.; Evall, J.; Whitesides, G. M.; Nuzzo, R. G. J. Am. Chem. Soc. 1989, 111, 321. (11) Poirier, G. E.; Tarlov, M. J.; Rushmeier H. E. Langmuir 1994, 10, 3383-86.
© 1997 American Chemical Society
Alkanethiol Adsorption on Au
due to dynamic constraints. Obviously, the observations of ref 6 implicitly carry information on the dynamics of gas phase adsorption. Other authors have tried to approach a better understanding of alkanethiol adsorption by theoretical means. Sellers et al.13 have provided a potential surface for the interaction of alkanethiols with gold surfaces in the reactively adsorbed thiolate state in order to make possible computer simulation studies of alkanethiolate films that are completed, i.e., after chemisorption and self-assembly have occurred. A second paper by Sellers is devoted to the investigation of the molecularly chemisorbed state of alkanethiols.14 This latter state is the precursor state of alkanethiol before it can reactively adsorb at the Au(111) substrate. In principle the knowledge of the precursor in addition to the adsorbed state should make possible a computer simulation not only of the completed film but of the adsorption process itself. However, as Sellers14 rightfully points out, the time scale of self-assembly is by far too large to be tractable by standard molecular dynamics techniques. In the present paper I make use of recently developed nonstandard molecular dynamics (MD) methods. They allow additional information on processes which take place on a time scale orders of magnitude greater than amenable in a direct way to molecular dynamics simulations. To visualize the situation, one has to realize that molecular dynamics studies are carried out on systems with a few thousand centers over periods of rarely more than a nanosecond, usually much shorter periods. On the other hand, as will be shown in the following for the adsorption of methanethiol on Au(111), one encounters elementary processes with a characteristic time of up to a few microseconds. The comparison between the simulation results and experimental data will even lead to the determination of characteristic times as long as 10-2 s. The method of “constrained dynamics” which has been devised during the last few years by several authors and concisely described in the literature,15-18 is suitable to solve the problem. The general idea is to define a reaction coordinate b rRC and to constrain the system with respect to this coordinate. When equilibrium under the constraint bRC) along the reaction is reached the mean force F BRC(r coordinate is sampled. With a speed low enough not to violate the equilibrium condition, the system is moved along the reaction coordinate. Integration of the force bRC) along the reaction coordinate allows computation F BRC(r of the work needed to push the system from the starting rt configuration b rs to the final configuration b
∫brbr 〈FBRC(rbRC)〉 drbRC
∆A(r bs,r bt) ) -
t
s
The reaction coordinate not only can be the length of a path along which a particle is to be transported,15 but can describe a rather complex conformational transformation of a macromolecule. In the latter context Schlitter has introduced the term “targeted dynamics”.16 In a complex (12) Alves, C. A.; Smith, E. L.; Porter M. D. J. Am. Chem. Soc. 1992, 114, 1222-1227. (13) Sellers, H.; Ulman, A.; Shnidman, Y.; Eilers, J. E. J. Am. Chem. Soc. 1993, 115, 9389-9401. (14) Sellers, H. Surf. Sci. 1993, 294, 99-107. (15) Mosell, T.; Schrimpf, G.; Hahn, C.; Brickmann, J. J. Phys. Chem. 1996, 100, 4571-81, 4582-90. (16) Schlitter, J.; Engels, M.; Kru¨ger, P.; Jacoby, E.; Wollmer, A. Mol. Simul. 1993, 10, 291-308. (17) Schlitter, J.; Steinhoff, H.-J. Parametrization of Conformational Transitions by the Distance in Configuration Space: Free Energy Profile. Preprint. (18) Mu¨lders, T.; Kru¨ger P.; Swegat, W.; Schlitter, J. J. Chem. Phys. 1996, 104, 4869-70.
Langmuir, Vol. 13, No. 15, 1997 3991
case like Schlitter’s the definition of the target conformation is unambiguous, but the choice of the reaction coordinate is not. Schlitter16 proposes a rather general way to define the reaction coordinate by a constraint as little perturbing as possible. The identification of the integral of mean force with the thermodynamic concept of free energy is not obvious, but generally accepted.15 Recently Schlitter has given a detailed justification.18 The question of which free energy is determined in this way, Gibbs energy G or Helmholtz energy A, simply depends on the side conditions. In the present case no external pressure is applied, which results in an identity between both quantities. In the above equation I have used the letter A (Helmholtz energy) in order to be in line with other authors.15,17 The paper is organized as follows. The following section describes the physical states a molecule encounters during the adsorption process. The rate constants for the transition between these states are the main results of this paper. Further, one has to define precisely which parts of the rather complex adsorption process can be treated presently in detail and which have to remain a black box. I would like to state at this point that rate constants of all processes can be derived. In the subsequent section information on involved interaction potentials is collected as found in the literature. A section follows with the results of the computer simulation calculations. The determination of rate constants as function of coverage and the construction of the full picture is given in the next section. An outlook on the perspectives, which the present paper on methanethiol opens to the treatment of higher homologues and finally to the treatment of adsorption from solution, completes this work. 2. States Relevant for Adsorption and Transitions between Them For the purpose of simulating the gas phase adsorption of methanethiol on Au(111) I distinguish four states: state 0, the methanethiol molecule is in the gas phase state 1, the methanethiol is in contact with already reactively adsorbed molecules, but its HS group is still separated from the substrate (physisorbed precursor state) state 2, the methanethiol molecule has its HS group in close contact with the surface, i.e., the HS group vibrates around the equilibrium position which is ∼3.2 Å in front of the surface (chemisorbed precursor state) state 3, the methanethiol is reactively adsorbed at the surface. The H-S bond is broken and the S atom vibrates around its equilibrium position which is ∼1.9 Å in front of the surface States 1 and 2 are distinguished only by the position of the SH group. At low coverage state 1 loses its meaning in that the approaching molecule cansand doessget immediately into contact with the substrate. However, a molecule approaching a precovered surface will at first have to accommodate on top of the organic film, i.e., into a typical physisorbed precursor state. Only then can the HS group penetrate the film and get into direct contact with the substrate which is typical for the molecularly chemisorbed molecule. States 0 through 2 are described by the same interaction potential with the substrate.14 We will see that for a short molecule like methanethiol state 1 is of relevance only for coverages above 0.9 monolayer (ML). In state 3 we do not deal anymore with the original species in that a chemical reaction must have taken place. It is assumed to be described by the net reaction
3992 Langmuir, Vol. 13, No. 15, 1997
metal + HSCH3 f metal-SCH3 + 1/2H2
Morgner
(1)
Whichever the details of the reactive adsorption process are, we describe the motion of the moleculesonce being in state 3sby the metal-SCH3 potential surface provided by Sellers.13 I will not deal explicitly with the chemical reaction (1). The reason is that the process is not necessarily unique, several reaction channels are discussed14 and observed indirectly by way of reaction products.19 In my model description of adsorption the transition between the chemisorbed state 2 and the reactively adsorbed state 3 is treated only by a rate constant k23. The chances of an adsorbing molecule to finally end in the reactively adsorbed state 3 (and thus to increase the coverage of the surface) depends on the period of time that it spends in the chemisorbed state 2. We will see later that the transition from state 2 to state 3 is the bottleneck for adsorption. The reactive sticking coefficient is mainly determined by the competition between the process 2 f 3 on one side and 2 f 1 or 2 f 0 on the other side. We will see later that only at high coverage the sticking coefficient is further influenced by the interplay between transition 1 f 0 (desorption out of physisorbed state) and transition 1 f 2. At low coverage the molecule passes through state 1 so swiftly that state 1 can be disregarded. In the following, the rate constant for transition from state i to state j is named kij. The main contribution of the present paper is to determine the rate constant k20 as function of coverage and at higher coverage the correction caused by the influence of state 1. Once this information is available, it is possible to determine the rate constant k23 by comparison with experiment. The reverse rate constant k32 is considered to be negligible at room temperature. The strategy for determining k23 runs via the sticking coefficient S. At low coverage state 1 is negligible. At high coverage the probability p that the adsorbing molecule reaches state 2 has to be considered. A molecule adsorbed into state 2 has two possibilities, reactive adsorption into state 3 or desorption into the gas phase. The probability of reactive adsorption is then
S)p
k23 k23 + k20
(2)
which is to be identified with the reactive sticking coefficient. k20 and p will be determined as a function of coverage. For θ e 0.9 the probability p turns out to be unity. The rate constant k23 is considered to be independent of coverage. This is done due to lack of better knowledge in the first place, but one argument is strongly in favor of this assumption. Reaction 1 leads primarily to chemisorbed atomic hydrogen which could block the adsorption of further molecules. However, it is known from energetic considerations6 as well as from direct observation that hydrogen desorbes very fast molecularly. We are left with the possibility that increasing coverage influences the rate of transition 2 f 3 dynamically. This cannot be excluded, but it is expected to be a small effect in case of reactive sticking via the “hydrogen channel” conceived here.14 k23 is not known from literature. Therefore one has to guess a trial value for k23. The working formula for the reactive sticking coefficient as function of coverage θ is then given as (19) Jaffey, D. M.; Madix, R. J. Surf. Sci. 1994, 311, 159-171.
S(θ) ) p(θ)
k23 k23 + k20(θ)
(3)
Now we are in the position to determine the relation between coverage θ and exposition Pt, where P denotes gas pressure and t time. The influx of molecules to the surface per unit area and per unit time at a given pressure P is given by
J ) CP
(4)
where
C)
NA (2π MRT)1/2
and NAis Avogadro’s constant, R is the gas constant, T is temperature, and M is the mass of 1 mol of molecule. With nML being the number of adsorbed molecules per unit area in a monolayer, we have
dθ nML ) S(θ) C d(Pt)
(5)
Integration gives
(Pt)(θ) )
nML C
dθ′ ∫0θ S(θ′)
(6)
This function gives the necessary exposition Pt in order to achieve the coverage θ. From experiment6 it is known that Pt > 104 langmuirs is needed to create a full monolayer of methanthiol on Au(111). Accordingly, one can adjust k23 to the value which reproduces the experimental result. In principle, it is possible to go one step further and determine experimentally the relation between exposure Pt and coverage θ. Then it should be possible to determine the rate constant k23 for any coverage and thus carry out a test on the assumption made here that k23 is independent of coverage. 3. Model Description of the System The methanethiol molecule is treated as being composed of the two pseudoatoms HS and CH3 coupled by a harmonic force. The rate constant for desorption out of state 2 (molecularly chemisorbed HSCH3) is computed for a rigid surface. The interaction of the methyl group CH3 depends always only on the distance z from the surface. If the molecule is in state 1 or 2, the same holds for the HS group. In the reactively adsorbed state the HS group is to be replaced by an S atom. Its interaction with the surface is carefully modeled according to ref 13 in order to describe not only the potential energy dependence on the distance z but also the lateral corrugation. The potential minimum is taken on at a hollow site on the Au(111) surface whereas the top site has a binding energy which is ∼0.260 eV smaller than in the hollow site. The data from ref 13 are reproduced in our model ansatz. A spline function allows for a smooth variation of potential energy within the unit cell. The parameters are given in the subsequent section. Taking into account the lateral corrugation for state 3 is important for the simulation of the self-ordering process on the surface. This is demonstrated by the comparison of alkanethiol adsorption on silver(111) and gold(111). The lateral order and film densities differ markedly even though the substrate lattices are almost identical and the interchain interaction is the same whereas the lateral corrugation is smaller in case of Ag(111). Therefore, it is worthwhile to put effort
Alkanethiol Adsorption on Au
into the modeling of state 3. In the course of film growth the adsorbing molecules accumulate in the reactively adsorbed state (state 3). The higher the coverage, the more the dynamical behavior of the already adsorbed molecules influences the fate of a newly adsorbed molecule which is still in state 1 or state 2. On the other hand, the lateral dependence of the HS surface interaction in state 1 or 2 has little influence. However, later on we will encounter one particular situation where the neglect of lateral corrugation in state 2 is necessarily unrealistic. This happens at very low coverage: it turns out that the desorption probability of a molecule in state 2 depends strongly upon whether it can bind to an already adsorbed molecule or not. Thus, the question arises whether the adsorbing molecule in state 2 has a chance to find a reactively adsorbed molecule before desorption. Since reactively adsorbed molecules at 300 K have a very low lateral diffusion constant, the encounter of the adsorbing molecule with an already sticked molecule depends only on the lateral mobility of the first. Thus, it is of interest to characterize the lateral motion of a molecule in state 2. This, however, cannot be done on a surface which is modeled as rigid and without corrugation. It is obvious that one has to introduce an additional scheme for the interaction of a molecularly chemisorbed methanthiol molecule with an Au(111) surface. One way to go would be the explicit incorporation of lateral corrugation as done for state 3. Another way could be to leave the model of a rigid surface and to simulate the gold surface by a cluster of suitably arranged gold atoms. Even though this represents increased computational effort, I have chosen the latter alternative. Several reasons are in favor of this choice: the nonrigid surface should be even more realistic in modeling the lateral mobility of the molecule whereas a corrugation is automatically introduced by a pair interaction potential between the pseudoatoms of the molecules and the substrate atoms. The choice of parameters for these pair potentials and the potential between the gold atoms is discussed below. Further, the representation of the substrate surface by a cluster of atoms (in this case 192) allows observation of the dissipation of the translational energy of the incoming molecule into the surface. Accordingly, one can learn from the molecule-cluster interaction whether the approaching molecule gets rid of its energy sufficiently fast to guarantee sticking or whether reflection occurs. The answer is that methanethiol is found to stick always and that it remains on the Au cluster for long times. On the other hand, it is obvious and has been found in calculations that within the model of rigid surface an approaching small molecule like methanethiol is reflected inevitably and almost instantly. Therefore, the possibility for energy dissipation (i.e., temperature control) must be explicitly built into the calculation when employing the rigid surface. The question is then how to choose the adequate relaxation time for temperature control. This number can easily be derived by comparison with the more realistic cluster model. The time step in all molecular dynamics calculations has been set to 2 fs. The speed along the reaction coordinate during “constrained dynamics” must be chosen sufficiently low in order not to violate the condition of bRC). After equilibrium when measuring the mean force F BRC(r some trial runs the motion along the reaction coordinate has been carried out with a velocitiy of 1.5 × 10-4 Å/fs. In the next section the analytical forms and their parameters are collected for the interaction of the molecules with the rigid surface and the atoms of the cluster and for the interaction between the molecules.
Langmuir, Vol. 13, No. 15, 1997 3993 Table 1. Parameters of LJ Pair Potentials (C6)1/2, (C12)1/2 (C6)1/2 (C12)1/2, species (eV Å12)1/2 (eV Å6)1/2 ((kJ/mol) Å12)1/2 ((kJ/mol) Å6)1/2 HSa CH3b a
417.17 618.08
9.32 10.38
4097.7 6071.2
91.547 101.96
Reference 20. b Reference 21.
Table 2. Harmonic Potential between S(HS) and CH3a force constant
a
distance, Å
eV/Å2
(kJ/mol-1)/Å2
1.817
18.85
0.195
Sellers et al.13 sp3 hybridization.
4. Potentials and Parameters The interaction between molecules is described by Lennard-Jones potentials the parameters of which are taken from ref 20 and 21 and comprised in Table 1. The energy unit used in the computer code is the electronvolt. For easier comparison with other work the parameters are given in SI units as well. The interaction between the pseudoatoms is treated as harmonic potential. In principle the force constant and the equilibrium distance should depend on the chemical environment of the molecule, i.e., the parameters should be different for the molecule in gas phase, in the molecularly adsorbed and reactively adsorbed state. However, since the bond order is the same in all states treated, the effect may not be overly severe. The parameters found in the literature (ref 13) for the reactively adsorbed state were employed throughout (Table 2). The interaction of the pseudoatoms with the surface is described by the analytical form found in ref 22
V(z) )
C12 (z - z0)
C3 12
-
(z - z0)3
(7)
Table 3 contains the parameters for expression 7 as well as the equilibrium distances ze and the potential well depths Ve derived from these parameters. The minimum energy position of methanthiol in the molecularly chemisorbed state (state 2) at the Au(111) surface has been calculated from the parameters in Table 3 by lowering the temperature in a MD run down to 1 K. The distances of the SH group and the CH3 groups are then zSH ) 3.197 Å and zCH3 ) 3.515 Å. The angle between the surface normal and the HS-CH3 axis turns out to be 100.1°. Thus, the chosen parameters reproduce well the minimum energy position from ref 14. The minimum energy is calculated as -0.5125 eV ) 49.5 kJ mol-1, which is 7% smaller than the theoretical value from ref 14, but very close to the experimental value quoted in ref 14. The entries for species S in Table 3 refer to the reactively chemisorbed state. The parameters are chosen so as to reproduce the ab initio data13 for this situation. As discussed in the previous section the sulfur-surface interaction potential is considered to depend on the lateral position. The strongest attraction is experienced if the sulfur atom sits on the hollow site. In the top site the binding energy is distinctly smaller. This is modeled by a C12 coefficient which is larger by a factor of 1.55 as against the C12 coefficient at the top site. The C3-coefficients are identical in both cases. The necessary interpolation between top and hollow site is done in the following way. (20) Van Gusteren, W. F.; Berendson, H. J. C. GROMOS library manual; Biomos: Groningen, 1987. (21) Bo¨cker, J.; Schlenkrich, M.; Bopp, P.; Brickmann, J. J. Phys. Chem. 1992, 96, 9915-22. (22) Hautman, J.; Klein, M. L. J. Chem. Phys. 1989, 91, 4994-5001.
3994 Langmuir, Vol. 13, No. 15, 1997
Morgner
Table 3. Interaction of HS, S, and CH3 with the Rigid Surface C12
Ve
state
z 0, Å
eV/Å12
(kJ/mol-1)/Å12
eV/Å3
(kJ/mol-1)/Å3
ze, Å
eV
kJ mol-1
CH3a HSb Sc
2 2 3 top 3 hollow
0.86 0.86 0.105
2938.5 3900 832.0
283522 376293 80276
1.792 7.5 15.56
172.9 723.6 1501
3.515 3.197 1.935
0.0718 0.4407 1.907
6.93 42.52 184.0
0.105
1382.5
133391
15.56
1501
2.03
1.647
158.9
Sc
c
C3
species
a Hautman and Klein, ref 22. b Parameters adjusted to reproduce the equilibrium position and potential energy of the entire molecule.14 Parameters adjusted to reproduce the data given in ref 13.
Table 4. Coefficients of Spline Polynomial in Unit Cell fspline(χ,η) ) g4(χ2 + η2) + g5χη + g6(χ3 + η3) + g7(χ2η + χη2) + g8(χ4 + η4) coefficient
value as fraction
g4 g5 g6 g7 g8
287649/38404 29565/38404 -447309/76808 94041/38404 -142155/9601
Let (ax,ay), (bx,by) be the lattice vectors of the Au(111) unit cell with (ax ) 2.89 Å, ay ) 0 Å) and (bx ) 1.445 Å, by ) 2.5028 Å). Any position on the surface can then be described by the coordinates χ, η as χ(ax,ay) + η(bx,by). The conversion of the Cartesian coordinates x,y into χ, η is carried out by
() ( ) () a b χ ) ax bx η y y
-1
x y
The coordinates χ, η have the advantage that they are well adapted to the symmetry of the unit cell. The positions of the top and hollow sites are defined as
hollow site: top sites:
(χ, η) ) (0,0) (χ, η) ) (-1/3,-1/3); (-1/3,2/3); (2/3,-1/3); (2/3,2/3)
Within χ,η ∈[-1/3,2/3] the incomplete polynomial of 4th order
fspline(χ, η) ) g4(χ2 + η2) + g5χη + g6(χ3 + η3) + g7(χ2η + χη2) + g8(χ4 + η4) is chosen as spline function. It takes on the value 0 at the hollow site and unity at the top sites. The values of the coefficients are presented in Table 4. It should be noted that the corrugation of the potential surface experienced by the sulfur in the reactively adsorbed state leads to a force component parallel to the surface. Thus, the sulfur atom carries out vibrations not only along the surface normal but also parallel to the surface. Further, the adsorbed molecule can switch to another site, i.e., to another gold atom. Within the potential used the potential energy barrier against this process is 0.26 eV ) 25.09 kJ as can be seen from the difference between the entries for S(top) and S(hollow) in Table 3. Since this energy barrier is large compared to kBT at room temperature, a spontaneous site hopping hardly occurs in the simulation of methanethiol. However, under the influence of an adsorbing molecule (state 2), site hopping of a reactively adsorbed molecule can indeed be observed. In the reactively adsorbed state no explicit interaction potential between CH3 and the surface is used. The separation of the CH3 group from the surface is controlled by a harmonic bending potential. The bending coordinate
is the angle between the S surface and the S-C directions. For the hollow site Sellers et al.13 found two equilibrium values of the bending angle at 180° (sp sulfur hybridization) and at 104° (sp3 sulfur hybridization), the latter being lower in energy and having a much stiffer force constant. In the top site, only the sp3 geometry persists as equilibrium situation.13 Since the motions of the reactively adsorbed molecules havesexcept for very low coveragesa strong influence on the dynamics of adsorbing molecules, it appears important to incorporate the bending potential into the MD studies. The parameters which are employed to model the metal-S-CH3 bending potential are given in Table 5. The situation at the hollow site has to be described by two harmonic potentials, one with the minimum at 104° and the other at 180°. The minimum energy at 180° lies about 0.018 eV above the value at 104°. In the intermediate range (115°, 130°) a polynomial of third order is used to connect both curves smoothly. The same analytical ansatz is used for the top site. However, since no second potential minimum exists at the top site,13 the parameters must be chosen in such a way that a second minimum is avoided. Firstly, the potential energy at 180° is set to a rather high value of 0.42 eV and the force constant is made negative. Both quantities are adjusted so as to effectively remove the second minimum. Between the hollow site and the top site the above mentioned two dimensional spline function (Table 4) is used to interpolate the force constant and the potential energy at 180°. As mentioned in the preceding section, some properties of the adsorbing methanethiol have to be calculated via the interaction not with a rigid surface but with a cluster of gold atoms. For this purpose, all relevant pair potentials have to be established. In order to find Lennard-Jones pair potentials the following approach has been used which can easily be generalized to any other metal. The pair potential between gold atoms can be derived from the bulk compressibility κ which is related to the energy U of N atoms as function of the volume V they can fill.23
U(V) )
V4
3/2
b6
b12 -
; V2
κ)
1 b12 x2 b65/2
(8)
The relation between the coefficients b12, b6 and the equilibrium volume V0 is
V0 )
( ) 2b12 b6
1/2
Here we are interested in the interaction potential between pairs of atoms which will be expressed as a LennardJones potential (23) Kittel, Ch. Einfu¨ hrung in die Festko¨ rperphysik; Oldenbourg: Mu¨nchen, 1983.
Alkanethiol Adsorption on Au
Langmuir, Vol. 13, No. 15, 1997 3995 Table 5. Bending Potential Metal-S-CH3
site hollow top
valid between
potential energy
force constant
equilibrium angle, deg
deg
deg
eV
kJ mol-1
eV/rad2
(kJ/mol) /rad2
104 180 104b 180
0 130 0 130
115 180 115 180
0 0.01778a 0 0.42c
0 1.716 0 40.52
2.01a 0.125a 2.01a -0.625c
193.9 12.06 193.9 -60.3
sp3 sp sp3
a Parameters adjusted to reproduce the data given in ref 12. b For simplicity of computer code instead of 102° as given in ref 12. c Adjusted to remove second minimum without changing analytical form.
Table 6. Pair Potential for Au-Au
V(R) ) quantity
symbol
compressibility number density well depth
κ N/V0 4 σ C12 C6 f
force constant
12
C12
C6 -
12
R
R
C12 )
R12
C6 -
R6
) 4
12
6
((Rσ ) - (Rσ ) )
value
unit
5.774 × 10-12 5.9 × 1028 102.65 2.644 1.198 × 107 3.507 × 104 209.4
m2/N m-3 kJ mol-1 Å kJ mol-1‚Å12 kJ mol-1‚Å6 kJ mol-1‚Å2
unit
1.064
eV
124187.8 363.5 2.17
eV Å12 eV Å6 eV Å-2
C12
(9)
1 P N5σ12 2 12
b6 ) P6N3σ6 In case of a face-centered cubic (fcc) crystal the factors P12, P6 have the values23 P12 ) 12.13188; and P6 ) 14.4539 and the relation between volume per atom V/N and next neighbor distance R is
R3 ) x2 V/N The calculated parameters including the force constant at the minimum energy position are given in Table 6. The well depth of the pair potential is 1.064 eV whereas the binding energy per atom in the equilibrium position is evaluated as 2.29 eV according to eq 8. This latter value is to be compared with the heat of formation of gold which is 3.816 eV.24 The two quantities which are reproduced by the parameters in Table 6 are the equilibrium separation of nearest neighbors and the force constant which enter via the number density and the compressibility into the treatment. Since the functional forms (8) or (9) have only two independent parameters, one is not surprised that a third physical quantity, i.e., the heat of formation cannot be described at the same time. A function with higher degree of flexibility would be needed. For the present purpose, namely, the treatment of a gold cluster at 300 K, the determined pair potential is well suited since the cluster will remain near its equilibrium configuration. However, it is obvious that the potential would lead to wrong results if employed at much higher temperatures, e.g., near the melting point. The remaining task is to find pair potentials between the pseudoatoms HS and CH3 and the Au atoms. This attempt must be guided by the requirement that the pair (24) Weast, R. C. Ed. CRC Handbook of Chemistry and Physics, 67 ed.; CRC Press, Inc.: Boca Raton, FL, 1986.
C6
pair
eV Å12
kJ mol-1 Å12
eV Å6
kJ mol-1 Å6
HS-Au CH3-Au
323206.8 207391.0
3.118 × 107 2.001 × 107
284.4 83.74
2.744 × 104 8.079 × 103
HS-Cluster of 192 Au Atoms
The relation between expressions 8 and 9 is given by
b12 )
value
Table 7. Pair Potentials HS-Au and CH3-Au
6
((Rσ ) - (Rσ ) )
W(R) ) 4
6
Ve site
ze,a Å
eV
kJ mol-1
hollow top
3.1 3.3
0.449 0.365
43.3 35.2
CH3-Cluster of 192 Au Atoms Ve site
ze,a Å
eV
kJ mol-1
hollow top
3.55 3.75
0.078 0.069
7.5 6.7
a Measured from the center of mass of the topmost layer of Au atoms.
potentials reproduce the already known interaction potentials between HS, CH3, and the rigid surface; cf. Table 3. As an analytical form, again a Lennard-Jones potential has been chosen. The parameters are listed in Table 7. Further, the effective interaction potentials with the Au(111) surface, as obtained by summation over all 192 Au-atoms of the gold cluster, are characterized by their equilibrium distances ze and their well depths Ve. A distinction is made between approaches towards hollow sites and top sites of the gold cluster. Since the interaction of the pseudoatoms with the rigid surface (Table 3, state 2) is modeled without lateral corrugation, one has to compare the means of top and hollow sites in Table 7 with the entries of Table 3. The agreement is reasonable. Sellers13 reports for HSCH3 a binding energy to the Au(111) surface which is higher at top sites than at hollow sites. Inspection of Table 7 shows that the LJ potential ansatz is not flexible enough to reproduce this result. The attraction is always stronger at the hollow site due to the rather hard repulsive wall of the LJ potential. For the present application this may not be overly severe because the purpose of the methanethiol Au-cluster simulation is to model energy dissipation and lateral diffusion which are controlled by the possibilities for energy/momentum exchange between molecule and surface and by the coupling of motion perpendicular and parallel to the
3996 Langmuir, Vol. 13, No. 15, 1997
Figure 1. Adsorption of methanethiol at an Au(111) surface modeled by a cluster of 192 Au atoms: dashed line, instantaneous temperature of the system, the averaged temperature is indicated; full line, z-position of the HS group of the adsorbing methanethiol molecule, the average z-position of the gold atoms of the first layer is at z ) 1.82 Å.
surface due to corrugation. Both mechanisms exist in the model used. With respect to the quantitative evaluation of lateral diffusion the amount of energy corrugation could be meaningful. The value found in the ab initio calculation13 is 0.143 eV whereas only 0.084 eV is obtained in the present model. Thus, one has to consider the possibility that the lateral mobility of methanethiol on the Au cluster is too large and therefore to be considered as upper bound to a realistic value. 5. Results Adsorption of HSCH3 on Au Cluster. The gold cluster consists of 192 atoms arranged in four layers of 48 atoms. Periodic boundary conditions are applied in the surface plane. The box lengths are 17.34 Å and 20.0224 Å in x and y directions. The Au atoms are arranged to simulate a Au(111) surface. The methanethiol molecule approaches along the z-direction. The translational energy of the methanethiol corresponds to 300 K at large distances. In addition it gains about 0.5 eV kinetic energy due to the attractive molecule-surface interaction. The goal of the simulation is to study whether dissipation of energy into the gold surface is fast enough in order to avoid direct reflection of the impinging molecule. The size of the gold cluster is so large that it can absorb the energy without undue temperature increase. The estimate amounts to ∆T e 15 K to be compared with the starting temperature of 300 K. In fact, a temperature rise of ∼10 K is observed within a time of about 1 ps during contact between the molecule and the cluster. The development of the temperature of the system over a period of 10 ps is shown in Figure 1. In order to visualize the relation to the adsorption process the z position of the HS group of the methanethiol molecule is shown as well. No temperature control is applied during this run. The simulation of the adsorption process has been repeated three times using different starting conditions of the molecule. The result was always similar. After adsorption the development of the system was followed for periods up to 800 ps. Spontaneous desorption was never observed. From the constrained dynamics studies described below this will be quite understandable: the lifetime against
Morgner
Figure 2. Lateral diffusion of methanethiol on (111) surface of gold cluster. The area probed by the molecule for the presence of reactively sticking molecules is shown as function of time after adsorption. The straight line represents the linear regression. Further explanation is given in the text.
desorption of an isolated methanethiol on Au(111) turns out to be much longer than the investigated period, namely, 0.6 µs. The important message from the above described simulation is that the methanethiol is not reflected from the gold surface but can deliver its (thermal) kinetic energy sufficiently fast to the substrate. This feature must then be built into the simulation if the adsorption onto a rigid model surface is to be studied. Lateral Mobility of HSCH3 on Au Cluster. For later use the question is investigated which area is explored by the adsorbed HSCH3 molecule in the course of lateral diffusion. As mentioned above, the goal is to determine the probability that the molecularly chemisorbed methanethiol finds an already reactively adsorbed SCH3. The word “finding” has to be specified: it is assumed that a lateral separation of 5 Å or less leads to bonding due to van der Waals attraction. Thus, the diffusing molecule explores at any time an area of πR2 with R ) 5 Å. The area Q(t) explored by the molecule in the course of time can be obtained if a circle of radius R is pushed along the trajectory (precisely the projection of the trajectory onto the surface plane). All points touched in this way contribute to Q. Since many overlaps occur during this procedure, an algorithm for computing Q from the trajectory would not be trivial. On the other hand, if understood as a graphical task, the solution is very simple: at any point of the trajectory a colored circle of radius R is drawn and the total number of colored points is then a measure for the desired quantity Q. Exactly this procedure has been applied by drawing the circles onto a computer screen and using the graphic tools of a suitable programming language for evaluation.25 The area Q(t) as a function of time is plotted in Figure 2. As expected, the data points can well be fitted by a straight line
Q(t) ) 447.75 + 5.4792t
(10)
where Q is in Å2 and t is in ps. The intersection with the ordinate is 447 Å2, which is much larger than the area πR2 ) 78.54 Å. Obviously, the temperature of the (25) Microsoft Visual Basic.
Alkanethiol Adsorption on Au
Langmuir, Vol. 13, No. 15, 1997 3997
adsorbing molecules needs a few picoseconds before the temperature is fully equilibrated with the substrate causing an initial “hot phase” with respect to lateral diffusion. Adsorption into State 2 at Rigid Substrate. A new molecule is created at distance z ) 25 Å above surface in its equilibrium conformation and the atoms are given velocity components according to the temperature of 300 K. Then the velocity components toward the surface (zdirection) are readjusted to give the center of mass of the molecule a velocity toward the surface (300 K). The molecule travels then for several picoseconds towards the surface and hits at a position which depends on the randomly chosen initial position and the randomly chosen velocity components in the x, y direction. When the first molecule arrives at the bare surface, it is reflected almost instantly if no temperature control is applied. This is easily understood: the surface is rigid (as opposed to the Au cluster) and the only degree of freedom into which the translational energy can flow is the vibration and rotation of the molecule itself. From the adsorption of methanethiol on a Au cluster we know that the dissipation of energy into the cluster (i.e., the surface) is fast enough to prohibit direct reflection. The rigid surface cannot do this, therefore energy dissipation must be explicitly built in. This is done by means of temperature control. The crucial quantity is the relaxation time. Several test calculations show that relaxation times down to 5 ps are still too long to prohibit reflection reliably. The final choice of the relaxation time was 2 ps. Essmann and Geiger26 used a relaxation time of 0.5 ps to achieve adsorption of H2O on a solid substrate. In the present case the longest relaxation time that prohibits direct reflection is used. Temperature control with the relaxation time of 2 ps is applied in all calculations at the rigid surface. This procedure is repeated until the number Nmax of molecules which corresponds to a x3×x3 monolayer is adsorbed. Calculations are carried out for two system sizes: ensemble I, Nmax ) 12, box length in x-direction 17.34 Å and y-direction 15.0168 Å; ensemble II, Nmax ) 48, box length in x-direction 34.63 Å and y-direction 30.0336 Å, the value of the nearest neighbor distance between Au atoms being 2.89Å. The cutoff radius for LJ interaction is set to rcut ) 8 Å. This is larger than half of the y box length in ensemble I. The comparison with ensemble II shows that this has no effect on the evaluated quantities; cf. the entries for the free energy of desorption in Table 8. The atom-surface interaction is cut off at zcut) 12 Å. Desorption from State 2 into Gas Phase. On the time scale of a MD simulation desorption out of the chemisorbed state (state 2) is not observed. Indeed, we will see that the lifetime against desorption out of state 2 is of the order of a microsecond which is orders of magnitude larger than the time covered by a standard MD calculation. It is obvious that a nonstandard technique must be used. We start with the expression
( )
kd ) Ad exp -
Edes kBT
(11)
for the desorption rate constant where Edes is the work needed to remove the molecule from the surface whereas the frequency prefactor Ad represents the average vibrational frequency of the system along the reaction coordinate which in our case is simply the distance zcm between the center of mass of the desorbing molecule and the surface. (26) Essmann, U.; Geiger, A. J. Chem. Phys. 1995, 103, 4678-92.
The key quantity is Edes. Techniques to compute it have recently been described in the literature.15-18 I follow the strategy to compute the free energy profile along the reaction coordinate. In the present case the definition of the reaction coordinate is easily possible: any valid starting configuration b rs is defined by the requirement that the desorbing molecule is in state 2 whereas all other molecules are firmly bound to the surface in state 3. State 2 is defined by the equilibrium distance of the HS group in front of the surface. The molecule is considered to be in state 2 if its HS group vibrates around the equilibrium distance ze(HS) ) 3.2 Å for a period as long as 40 ps. In some runs this time interval has been reduced to 4 ps with no noticeable effect. The final or target configuration b rt is defined by the observed molecule being irrevocably separated from the surface. The condition that the center of mass of the molecule has reached a distance zcm ) zt ) 15 Å from the surface proved sufficient because the force between the molecule and the surface has dropped to zero. The reaction coordinate itself is then defined as zcm, the distance between the center of mass of the molecule and the surface. The starting value of the reaction coordinate zs is consequently the value of zcm at the instant when the procedure of targeted/constrained dynamics is begun. Note that this value may vary from run to run for two reasons: first, the definition of the starting configuration is not directly given via zcm but via z(HS), i.e., the z-position of the HS group of the molecule, and second, the state is not defined by a fixed value of z(HS) but HS vibrating around its equilibrium position. The free energy profile
∆A(zcm) ) -
∫zz s
〈FRC(z)〉 dz
cm
(12)
is a function that essentially increases monotonically with zcm. Only a few small wiggles are superimposed to this general behavior. At high coverage, one reproducibly finds a negative slope at a distance from the surface which corresponds to the outer boundary of the reactively adsorbed molecules. Apparently, this signals the increase in entropy for the system when the molecule is pulled out of the molecular film. We must state, however, that the true gain in entropy upon desorption is much larger than that which corresponds to this small feature in the free energy profile. Apparently the gain in entropy occurs along the major part of the reaction pathway and is masked by the increase in energy. We observe that at large distances the free energy profile reaches its asymptotic value rather smoothly. No intermediate energy barrier is encountered. This indicates that a localized transition state does not exist. This is fortunate since it avoids the problem to determine the “transmission coefficient”. It is obvious that the desorbing molecule whenever it reaches the asymptotic range with zcm > zt must have a velocity component of its center of mass which is directed away from the surface. Since at zcm > zt all forces between surface and molecule are zero the molecule continues its trajectory to infinity. The asymptotic value of the free energy is employed as desorption energy
Edes ) ∆A(zcm ) zt)
(13)
Sincesas outlined abovesthe definition of the starting configuration allows for a lot of different configurations, we repeat the procedure at least 10 times with significantly different starting configurations. We consider a new starting configuration to be significantly different if it has evolved in unconstrained dynamics for at least 4 ps, in most runs even for 40 ps, or if a completely new ensemble has been constructed. Up to 30 values for Edes were
3998 Langmuir, Vol. 13, No. 15, 1997
Morgner
Table 8. Free Energy for Desorption of HSCH3 Molecules at 300 Ka no. of molecules
Edes mean, eV
std dev, eV
ensemble I: 12 SCH3 molecules per monolayer; box lengths of periodic boundary conditions (x,y) ) (17.34 Å, 15.0168 Å) 1 0.36 0.033 2 0.46 0.032 3 0.46 0.036 4 0.44 0.052 5 0.46 0.027 6 0.46 0.036 7 0.43 0.053 8 0.43 0.046 9 0.40 0.061 10 0.35 0.054 11 0.36 0.057 12 0.36 0.041 ensemble II: 48 SCH3 molecules per monolayer; box lengths of periodic boundary conditions (x,y) ) (34.68 Å, 30.0336 Å) 1 0.36 0.025 2 0.47 0.043 3 0.46 0.026 4 0.48 0.021 5 0.47 0.025 12 0.45 0.022 13 0.44 0.022 14 0.46 0.031 25 0.43 0.028 34 0.38 0.021 35 0.39 0.039 36 0.40 0.026 37 0.35 0.028 48 0.36 0.034 a
Figure 3. Lateral diffusion of HSCH3 on top of a dense methanethiolate layer. The area probed by the molecule for the presence of an empty site within the layer is shown as a function of time. The straight line represents the linear regression. Further explanation is given in the text.
Rate constant is given by kd ) Ad exp(-Edes/kBT).
collected in this way. The values given in Table 8 are the mean values 〈Edes〉 and the standard deviations to give an insight into the range of values encountered. To estimate the desorption rate constant we still need the frequency prefactor Ad. The prefactor has to reflect the vibrational motion of the observed molecule in the starting configuration along the reaction coordinate. Consequently we have evaluated the z-motion of the center of mass of the observed molecule in state 2 under unconstrained dynamics. By fast Fourier transform the frequency spectrum of zcm is computed.27 The mean value of this frequency spectrum is then employed as a value for Ad. This procedure is carried out for all periods of unconstrained motion that the system is allowed to perform before Edes is calculated by constrained dynamics. Accordingly, we have as many values for Ad as we have values for the desorption energy. No correlation between Ad and Edes was found. The mean values and standard deviations of Ad have been calculated for the smaller ensemble I (see Figure 4). Since the frequency factors Ad depend only very mildly on coverage, i.e. the number of the molecule considered, the procedure has not been applied to the larger ensemble II. Lateral Mobility within and Desorption from State 1. At high coverage the state 1 begins to play a role. The adsorbing molecule may have to remain within this physisorbed state for some time before it finds an empty site between the already reactively adsorbed -SCH3 molecules and reaches state 2. Its lateral mobility in state 1 as well as the lifetime against desorption are of interest. The computation of these quantities is performed by using ensemble I with all 12 molecules being adsorbed (27) Subroutine ‘realft’ from: Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical recipes in Fortran. The Art of Scientific Computing 2nd ed.; Cambridge University Press: Cambridge, 1992.
Figure 4. Frequency of the center of mass of methanethiol in state 2 along the surface normal as a function of the number N of the molecule. The calculation is carried out within ensemble I, i.e., 12 adsorbed molecules correspond to one monolayer. The linear regression yields Ad(N) ) ((2.14 × 1012-2.27 × 1010)N) s-1 for 1 e N e 12.
into state 3, i.e. a completed monolayer of methanethiolate is simulated. On the top of this ensemble an additional HSCH3 molecule is placed. Since the monolayer is complete, this molecule should not be able to penetrate into the layer of adsorbed molecules. The only way to leave state 1 should be desorption into the gas phase. The molecule in state 1 has been studied by unconstrained molecular dynamics simulation over a time interval of altogether 800 ps. During this time, penetration into the existing monolayers never occurred but spontaneous desorption was observed several times. Evaluation of the data yields a rate constant for desorption of k10 ) 1.358 × 1010s-1, which is equivalent to a lifetime against desorption τ10 ) 1/k10 ) 73.6 ps. Within one of the time periods which the molecule stayed in state 1 the trajectory of the HS group of the molecule was stored. In the same way as described above for HSCH3 at the cluster surface the area Q h (t) explored by the HS group could be evaluated. The development of Q h (t) over 27 ps is shown in Figure 3. The linear regression yields
Alkanethiol Adsorption on Au
Q h (t) ) (78.2 + 8.58t) Å2,
Langmuir, Vol. 13, No. 15, 1997 3999
t in ps
(14)
The value of Q h (t)0) is close to πR2 ) 78.5 Å2 based on the estimated range of interaction of R ) 5 Å. In contrast to adsorption onto the substrate itself, i.e., into the chemisorbed state 2, the adsorption into the physisorbed state 1 does not lead to an initial “hot phase”. 6. Discussion Rate Constant of Desorption as Function of Coverage. For comparison with experiment one would like to compute all quantities as a function of the coverage. Strictly speaking, the link between the number N of the adsorbing molecule and the coverage is yet to be established. Necessarily, the computer simulation relates to a finite (even though periodically repeated) area on which one monolayer corresponds to a maximum number of molecules Nmax. Thus, during adsorption of the Nth molecule the coverage changes from (N - 1)/Nmax to N/Nmax. Accordingly, any quantity computed in the simulation corresponds to a coverage θ which is (N - 1)/Nmax e θ e N/Nmax. In some special cases the choice is obvious: During adsorption of the first molecule one must set θ ) (N - 1)/Nmax. As long as the area of the periodic box is large enough to prevent interaction of the adsorbing molecule with itself, the first molecule represents the situation of a single molecule that adsorbs on a bare surface, i.e., θ ) 0. For the last molecule adsorbing, i.e., N ) Nmax, one has to choose θ ) N/Nmax ) 1 because this is adsorption into an isolated empty site which occurs at θ )1. For the intermediate cases a reasonable guess would be θ ) (N - 1/2)Nmax. This choice is employed in Figure 5 where the desorption energy Edes of methanethiol molecules from the chemisorbed state 2 is displayed as a function of coverage. One observes two plateaus at low and at high coverage and a distinct fall off around θ ) 0.6. An interesting observation is made at very low coverage. The extrapolation of the values down to θ ) 0 would yield a value of Edes ≈ 0.47 eV, but the value determined directly for N ) 1 is significantly lower for both ensembles. Such a discrepancy is unphysical, we expect that Edes (and thus the desorption rate constant and consequently the sticking coefficient) varies smoothly with θ. The way to treat smaller values of θ would be to treat ensembles with larger Nmax and to investigate whether Edes(N ) 1) and Edes(N ) 2) begin to deviate less for larger Nmax. This, however, would not be a very elegant method and further, we have no intrinsic criterion to tell which order of magnitude Nmax ought to be. Therefore the following strategy is followed: the different behavior between N ) 1 and N ) 2 appears to be due to the possibility for molecule N ) 2 to associate with the first reactively adsorbed molecule and thus be in a situation which differs from the first molecule. The question would then be: at which coverage does an adsorbing molecule have a good chance to encounter an already adsorbed molecule? This depends on two properties: the length of the time τdes during which the molecule remains in state 2 before it desorbes and the size of the area Q(τdes) that the molecule explores within τdes. The probability that the molecule does not find an adsorbed molecule is then
exp(-Q(τdes)θnML) where θ is the coverage and nMLis the number of adsorbed molecules per unit area in 1 monolayer.
Figure 5. Desorption of methanethiol out of state 2. The free energy of desorption is shown as a function of coverage. The dashed line is drawn to guide the eye.
Let kneigh be the desorption rate constant for a molecule that has found a neighbor and kisol be the desorption rate constant for a molecule without a neighbor. We obtain as effective rate constant for desorption
k(θ) ) exp(-Q(τdes)θnML)kisol + (1 - exp(-Q(τdes)θnML))kneigh (16) The relevant lifetime against desorption refers to a still isolated molecule, thus τdes ) 1/kisol and kisol is to be computed from Table 8 with the entries for N ) 1. One gets kisol ) Ad(N)1) exp(-Edes(N)1)/kBT) ) 1.72 × 106 s-1 and τdes ) 5.8 × 10-7 s. An estimate for the explored area Q as function of time has been derived from the methanethiol/Au cluster simulation as Q(t) ) (5.48t + 447) × 1012 Å2, with time t in seconds. It follows that the area that the adsorbing molecule can probe for the presence of an already sticked molecule is rather large, namely, Q(τdes) ) 3.2 × 106 Å2. With nML ) 0.0461 Å-2 one finds that the probability of encountering a sticked molecule becomes 0.5 at a very small coverage of θ ≈ 5 × 10-6. In section 4 it is discussed that the energy corrugation of the Au cluster surface is smaller by 0.059 eV than that obtained in Seller’s ab initio calculation. The Boltzmann factor would therefore be smaller by the factor exp(-0.059 eV/ kBT) ) 0.102. In consequence it could well be that the determined value for Q(t) is too large by an order of magnitude. If so, the point of switching over from the behavior of an isolated molecule to a molecule associated to a sticked molecule would shift from θ ≈ 5 × 10-6 to θ ≈ 5 × 10-5. We see that in any case the initial stage of adsorption which is characterized by the behavior of an isolated molecule in front of a bare surface is almost insignificant within the growth of a complete monolayer. The Reactive Sticking Coefficient. In order to determine the sticking coefficient at high coverage we have to determine the probability p that a molecule adsorbed into state 1 actually finds an empty site where it can proceed into state 2; cf. definition of p in eq 2. The probability p is given as
p ) 1 - exp(-Q h (τ10) (1 - θ) nML)
(17)
using the lifetime against desorption τ10, the explored area Q(τ10) within this period of time and the density of a monolayer nML.
4000 Langmuir, Vol. 13, No. 15, 1997
Morgner
Figure 6. Sticking coefficient and exposure as a function of coverage: full squares, sticking coefficient; open squares, exposure. Based on a rate constant for the transition from state 2 to state 3 of k23 ) 65 s-1.
The sticking coefficient S can then be evaluated for the whole range of coverages according to eq 3. As discussed in the second section the rate constant k23 can be determined via eqs 3 and 6 from the experimental information that adsorption of a monolayer requires an exposition of Pt > 104 langmuirs. Assuming that the exact value is Pt ) 1.5 × 104 langmuirs, the rate constant for transition from the chemisorbed state 2 into the reactively adsorbed state 3 is k23 ) 65 s-1. The sticking coefficient S(θ) computed on the basis of this value is shown in Figure 6. In addition, the exposure necessary to obtain a given coverage θ is shown as well. The slope of this latter function varies strongly between θ ) 0.7 and 0.8. A coverage of θ ) 0.7 can be effected already by a comparatively small exposure of 103 langmuirs. Adsorption beyond that point in order to complete the monolayer requires an additional exposure which is more than 14 times larger. Note, that this result is based on the assumption that the rate constant k23 is independent of coverage. Once experimental information on the relation between exposure and coverage is available, the validity of this assumption could be tested. The above quoted value of k23 ) 65 s-1 relies on the arbitrary decision that Pt > 104 langmuirs means 1.5 × 104 langmuirs. To get a better feeling for how crucial this decision is, Figure 7 shows how k23 would depend on a different choice for the monolayer exposure. Obviously, a precise measurement of the coverage as function of exposure would allow a rather precise estimate of k23 for all coverages. Within the framework of transition state theory, the rate constant k23 can be expressed as
k23 ) A κ exp(-Eact/kBT)
(18)
with the frequency prefactor A, transmission coefficient k, and the activation energy Eact. The route of the net reaction from state 2 to state 3 is not unambiguously known, but the characteristic values of Eact for different reaction paths have been discussed in the literature.14 For the “hydrogen channel” (see eq 1) Sellers quotes the theoretical value 70.7 kJ/mol and the experimental value 75.3 kJ/mol. We will try to obtain an estimate for the activation energy from the rate constant k23. For this purpose we have to guess a value for the frequency factor
Figure 7. Relation between the exposure necessary to built up a monolayer and the rate constant k23 for transition from the chemisorbed precrusor state to the reactively adsorbed state.
A. From the stretching vibration typical for hydrogen in organic matter, one estimates A ) 8.9 × 1013 s. For lack of better knowledge κ is replaced by unity. This means that the determined activation energy Eact ) 69.6 ( 1.2 kJ/mol is to be considered as the upper limit to the true value. The uncertainty in Eact is due to the uncertainty in k23. The apparent reasonable agreement may be taken as an indication that the present simulation widely reflects the real processes during adsorption. Entropy Gain during Desorption from the Chemisorbed State 2. The decrease of the free energy of desorption at higher coverage (cf. Figure 5) raises the question about the role of the entropy change ∆S in the process. From the equation
∆A ) ∆U - ∆S T
(19)
it follows that in addition to ∆A one has to compute the change of total internal energy ∆U during desorption in order to determine ∆S. Since the total energy U is the sum of kinetic and potential energies U ) Epot + Ekin, one has to evaluate
∆U ) Epot(after desorption) Epot(before desorption) + Ekin(after desorption) Ekin(before desorption) (20) As long as the temperature is kept constant, the kinetic energy should not change during desorption. Therefore, an equivalent expression should be
∆U ) Epot(after desorption) Epot(before desorption) (21) Both quantities, total energy and potential energy are available during any run. In practice, however, eq 21 is preferable for the following reason. For eq 20 to be correct, the quantities before desorption as well as after desorption had to be computed by unconstrained dynamics. As discussed in the previous section the evolution of the system is followed before desorption without constraints for a period of 40 ps. This is ample time to determine averaged values for potential and kinetic energy. The process of desorption is investigated by constrained dynamics and is considered complete as soon as the
Alkanethiol Adsorption on Au
Langmuir, Vol. 13, No. 15, 1997 4001
I. At higher coverage, i.e., for N g 10, the desorption entropy turns out to be rather close to the entropy of vaporization of ordinary liquids which is 85 J mol-1 K-1 according to Trouton’s rule.9 This may appear surprising at first glance. However, one has to keep in mind that the process studied here is not the desorption of a reactively adsorbed thiolate which is part of the self-assembled layer but of a nonreactively adsorbed chemisorbed species. One may conclude that surrounding thiolate molecules provide an environment to the HSCH3 molecule in the chemisorbed precursor state which is dynamically similar to an embedding liquid. 7. Summary and Outlook
Figure 8. Desorption of methanethiol out of the chemisorbed state. The calculation is carried out within ensemble I, i.e., 12 adsorbed molecules correspond to one monolayer. The number of molecules N is related to the coverage via (N - 1)/Nmax e θ e N/Nmax. Table 9. Entropy of Desorption for HSCH3 Molecules at 300 K ensemble I: 12 SCH3 molecules per monolayer; box lengths of periodic boundary conditions (x,y) ) (17.34 Å, 15.0168 Å) ∆U/eV no. of molecule mean std dev 1 2 3 4 5 6 7 8 9 10 11 12
0.45 0.45 0.51 0.48 0.51 0.53 0.53 0.54 0.55 0.60 0.62 0.60
0.007 0.030 0.066 0.054 0.061 0.064 0.063 0.086 0.090 0.090 0.118 0.066
∆S [eV/K] mean
std dev
0.00027 -0.00002 0.00019 0.00016 0.00017 0.00026 0.00034 0.00038 0.00052 0.00083 0.00085 0.00079
0.00010 0.00014 0.00023 0.00027 0.00021 0.00023 0.00030 0.00031 0.00037 0.00040 0.00047 0.00028
∆S [J/(mol K)] mean std dev 26 -2 18 15 16 25 33 36 50 80 82 76
10 14 23 26 20 22 29 30 36 39 46 27
desorbing molecule is outside the influence of the surface. The change of the free energy ∆A is computed during this part of the simulation. Obviously, it would require now an additional run without constraints to compute Epot and Ekin after desorption. This would consume additional computer time. Therefore, I have tried to make use of quantities already known. The question is whether the energies calculated under constraints during the late stage of desorption allow an estimate of the outcome of an unconstrained computation. Test calculations provide the following guideline: The larger the system the smaller is the influence of the constraint; this is as expected since the constraint concerns only one degree of freedom, irrespective of the number of molecules. The potential energy appears hardly affected whereas the kinetic energy reacts on the constraint. As a compromise, I have used eq 21 to determine the change of internal energy during desorption rather than performing new MD runs. The result should be very reliable at high coverage and somewhat less so at low coverage. The computed values for ∆U and for the entropy of desorption ∆S are given in Table 9. In Figure 8 the gain in entropy during desorption from the chemisorbed state 2 is shown. The calculations are carried out for ensemble
The process of gas phase adsorption of methanethiol on Au(111) is investigated by means of a molecular dynamics calculation, the working temperature being 300 K. Some of the elementary steps take orders of magnitude more time than can be handled by the standard MD technique. Therefore the tool of “constrained or targeted dynamics” has been used. The physical model invokes four different states in which an adsorbing molecule can be. The final state is that of reactively adsorbed thiolate which is named state 3. State 2 refers to the methanethiol molecule chemisorbed to the substrate via the HS group. State 0 defines the gas phase and state 1 the methanethiol physisorbed on top of an almost complete thiolate layer. It is obvious that state 1 does not play a role for coverages much below one monolayer. However, at coverages θ g 0.9 an approaching molecule first adsorbs into state 1. The probability for the molecule to get into state 2 depends on two quantities: the lifetime τ10 against desorption and the probability to find within this time an empty site where it can slip into state 2. For lower coverages the approaching molecule lands directly in state 2. There it experiences the competition between desorption back into the gas phase and the reactive adsorption into state 3. Accordingly, the desorption rate constant k20 and the reaction rate constant k23 must be determined. The desorption lifetime τ20 ) 1/k20 amounts to a few microseconds and, thus, requires the above mentioned special technique for its computation. Once this task is done for the full range of coverages, it is possible to determine the reactive rate constant k23 by comparison to available experimental data, namely, the exposure needed to form a dense and ordered monolayer of methanethiolates on Au(111). From the value of k23 it is possible to estimate the size of the corresponding activation energy. It agrees well with values derived from experiment and from ab initio data quoted in the literature. A key feature of the present ansatz is that the reaction 2 f 3 involves only the HS group of the molecule. The rest of the molecule has influence only as much as it influences the period of time the HS group can spend in state 2, which in turn controls the probability of transition into the final thiolate state 3. This has the important consequence that the reaction rate constant k23 is valid not only for methanethiol but also for higher homologues HS(CH2)mCH3 as well. The action of the HS group, of the methylene groups (CH2) and of the methyl group can be assessed by computer simulation techniques. Thus, the behavior in state 2 can be modeled for all alkanethiols and it can be calculated explicitly how the chain length influences the period of time the HS group can spend in the position typical for state 2. In consequence, the sticking coefficient as a function of coverage can be computed for all alkanethiols without any further information. Preliminary calculations predict monolayer exposures of only a few langmuirs for hexanethiol whereas the required exposure for nonanethiol appears to be larger
4002 Langmuir, Vol. 13, No. 15, 1997
by orders of magnitude. For comparison, the monolayer exposures reported in the literature are
