THE JOURNAL OF CHEMICAL PHYSICS 139, 174303 (2013)
Time domain simulations of chemical bonding effects in surface-enhanced spectroscopy Patrick Z. El-Khoury,a) Eric J. Bylaska, and Wayne P. Hessb) Physical Sciences Division, Pacific Northwest National Laboratory, P.O. Box 999, Richland, Washington 99352, USA
(Received 12 July 2013; accepted 16 October 2013; published online 4 November 2013) The atom-centered density-matrix propagation method is used to illustrate how time-dependent conformational changes affect the electronic structure and derived spectroscopic properties of a prototypical finite metal cluster-bound π -conjugated organic complex, Ag7 -benzenethiol. We establish that there is considerable conformational flexibility to the model structure, even at relatively low temperatures, which influences the predicted spectroscopic properties. Namely, the computed electron densities, dipoles, and polarizabilities are all dictated by torsional motion which controls the coupling between the π -framework of the chemisorbed molecular system and the cluster. © 2013 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4827455] I. INTRODUCTION
To account for the underlying dynamics in surfaceenhanced spectroscopy, it is essential to consider the distinct local environments in which different molecules reside. Electric field localization at defect sites in metal substrates induces a selective enhancement of the optical response of nearby molecules. This is a direct consequence of the inherent heterogeneity of typical metallic terrains over which molecular reporters are physi- or chemi-sorbed. In this regard, early works in the area of surface-enhanced spectroscopy1–4 inspired several ensuing studies,5–7 in which engineered hotspots were found to afford single molecule detection sensitivity. Understanding molecule-metal coupling in hybrid molecularmetallic constructs is a real challenge, both from the theoretical and experimental standpoints. It is a problem which lies at the interface of classical and quantum theories — each of which is individually tailored to solely describe either molecule or metal. A practical approach to rationalize chemical bonding effects in SERS is to simulate the local environment of a molecule using static tools of density functional theory (DFT) and its time-dependent variant (TD-DFT).8 In a pioneering work by Schatz and co-workers, the localization of electromagnetic fields at hotspots within a metallic substrate is identified as one of three cooperative mechanisms in surface-enhanced Raman spectroscopy (SERS).9 The other two are chemical in nature, namely (i) changes in the static polarizability of a molecule upon adsorption, and (ii) resonant enhancement due to metal-molecule charge transfer transitions. Of direct relevance to our present report is another early study by Aspuru-Guzik and co-workers, where the effects of chemical bonding on Raman scattering from benzenethiol (BT) chemisorbed on silver clusters was investigated using DFT calculations.10 Using different finite silver
a) [email protected] b) [email protected]
0021-9606/2013/139(17)/174303/5/$30.00
cluster models (Agn , n = 6–11), their work illustrated the need to consider various conformations of BT with respect to the finite cluster to correlate the metal-induced changes in the electronic structure of a molecule with SERS enhancement. Of the different effects governing normal (off-resonance) SERS enhancement, three were highlighted, namely (i) the relative orientation of the π framework with respect to the cluster, (ii) the local symmetry of the benzene ring, and (iii) the proximity of the different vibrational modes with respect to the binding sites. Our current work is driven by these observations and general interest in SERS. Rather than considering different clusters to emulate the aforementioned effects using static tools of quantum chemical simulations, we perform ab initio molecular dynamics (AIMD) simulations at the DFT level to systematically explore operative mechanisms in hotspot spectroscopy using one model system. Beyond the static picture, i.e., computing the optical properties starting at energy-minimized structures, the effect of time-dependent conformational and structural rearrangement in finite metal cluster-bound systems on their spectroscopic properties remains uninvestigated. Recognizing that the latter can at the very least be driven by thermal population under ambient laboratory conditions, we undertake the challenge of correlating the temporal evolution of molecular structure on a picosecond timescale with the derived electronic and vibrational structures in a prototypical model system, namely Ag7 -bound BT (PhS-Ag7 ). To this end, we perform AIMD simulations using the atom-centered density-matrix propagation (ADMP) method of Schlegel et al.11–13 The PhS-Ag7 system was chosen as a prototypical system because it has been well-characterized in previous theoretical10 and experimental works.14 This work extends previous theoretical studies by correlating time-dependent structural evolution with spectroscopic properties. We establish that conformational flexibility to our model structure dictates the computed electron densities, orbital energies (HOMO-LUMO gaps), dipole moments, and polarizabilities.
139, 174303-1
© 2013 AIP Publishing LLC
174303-2
El-Khoury, Bylaska, and Hess
J. Chem. Phys. 139, 174303 (2013)
II. COMPUTATIONAL METHODS
The B3LYP exchange-correlation functional and the LANL2DZ basis set are used throughout this study. Of the two possible isomers of PhS-Ag7 , we select the higher energy conformer (the vertex bound PhS-Ag7 structure shown in Figure 1) as a starting point for geometry optimization and frequency analysis. In principle, our calculations would allow for the interconversion between the two isomer forms where enough initial kinetic energy is partitioned into the isomerization coordinate. The harmonic IR and Raman spectra are represented as sums of Lorentzian functions that were empirically broadened by 15 cm−1 . In this concept report, the reduced basis set (266 basis functions) was chosen to decrease the computational burden of our rather resourceexhaustive simulations. However, in comparing our computed B3LYP/LANL2DZ properties at the minimum to their counterparts at the PBE0/def2-TZVP level of theory,10 it appears that the reduced representation is reasonable for the purpose of our current report, see Table I. Further validation of the level of theory used in the ADMP scheme follows. In Table I, we estimate the binding energy (Eb ) between BT and the Ag7 cluster is according to Eb (P hS − Ag7 ) = E(Ag7 ) + E(P hSH ) − E(P hS − Ag7 ) − 0.5E(H2 ).
(1)
We neglect basis set superposition errors and find that the computed Eb is reasonable, falling within the range of binding energies previously reported for BT bound to different sizes of silver clusters.10 The GAUSSIAN 0915 implementation of the ADMP methodology was used to carry out the AIMD simulations. The method is well-described in the original publications.11–13 Briefly, these are Car-Parrinello type AIMD simulations, carried out by propagating the density matrix and performed using atom centered Gaussians. The density matrices were obtained at each time step (t = 0.2 fs) from fully converged
TABLE I. A comparison of (i) selected structural parameters (bond lengths, angles in Ångstroms, degrees), (ii) binding energies (kcal mol−1 ), (iii) isotropic polarizabilities (a.u.), and (iv) HOMO-LUMO gaps (eV) in the PhS-Ag7 model computed in this work with their counterparts from Ref. 10. C2 -Ag denoted the shortest distance between the aromatic ring (a carbon atom adjacent to the C1 -S bond) and the cluster. The binding energy is calculated according to Eq. (1) in the main text. Source
Ag-S
C-S
This work 2.528 1.844 Reference 10–I 2.421 1.757 Reference 10–II 2.521/ 1.768 2.589
< C-S-Ag C2 -Ag 99.795 96.606 103.673
2.719 2.669 2.629
Eb
α iso
Egap
22.05 371.5 2.34 17.71 365.2 2.80 20.67 374.9 2.54
SCF calculations. The relatively small time step is chosen such that the total energy was converged to within ∼10−4 hartree. In these calculations, we use the same level of theory used for geometry optimization and harmonic frequency calculations, and employ the Cholesky decomposition scheme.16 We began by computing ADMP trajectories starting from the minimized structure. The initial nuclear velocities were randomly partitioned so that the total nuclear kinetic energies were 0.01, 0.05, 0.1, 0.15, and 0.2 hartree for five different trajectories of over 1 ps in total time duration. Next, to study the effects of large amplitude motion on the derived electronic and vibrational structures in the system, we identify a normal mode (ν 9 = 70 cm−1 ) of the vertex-bound structure that is best described as torsional motion along the C2 -C1 -S-Ag dihedral, see Figure 1(a). We translate the optimized structure by different increments along ν 9 , and use the resulting geometries as starting structures for another set of ADMP trajectories. Of the ADMP simulations performed, we analyze two representative trajectories in Figures 1 and 2. The first is initiated from the vertex-bound structure (henceforth designated as the reference trajectory). The starting structure for the second trajectory was obtained by displacing the same structure along the C2 -C1 -S-Ag dihedral bend, such that
0.6
(b)
Normalized Intensities (a.u.)
Relative energy (kcal mol-1)
(a) 0.5
0.4
0.3
0.2
0.1
0.0 0.0
0.1
0.2
0.3
Displacement along 9 (arbitary units)
0.4
0
500
1000 1500 2000 2500 3000 3500 Frequency (cm-1)
FIG. 1. (a) Change in energy for translating the minimized structure (0,0) by different increments along the C2 -C1 -S-Ag dihedral bend (torsion), a vector representation of which is shown in the inset. (b) Comparison of the Infrared and Raman spectra computed at the minimum to the Fourier transform of the nuclear velocity autocorrelation functions derived from ADMP simulations. The black trace (bottom) corresponds to the vibrational spectrum derived from the reference ADMP trajectory. Five frequency traces obtained by translating the minimum energy structure by different increments along the C2 -C1 -S-Ag dihedral bend at t = 0 were averaged to obtain the spectrum shown (dark yellow trace, second from the bottom).
J. Chem. Phys. 139, 174303 (2013)
Energy (Hartree)
-1262.26 TE - ref
TE - torsion
-1262.28 -1262.30 -1262.32 -1262.34
PE - ref
PE - torsion
(a)
-1262.36 0
200
400 600 time (fs)
800
C-S-Ag Angle (Degrees)
120 110 100 90 80
(c)
70 0
200
400
600
800
1.0
1000
time (fs)
(b)
0.8
FFT - ref FFT - torsion
0.6 0.4 0.2 0.0
1000
130
ref torsion
Normalized Amplitude (a.u.)
El-Khoury, Bylaska, and Hess
C-C-S-Ag Dihedral (Degrees)
174303-3
0
300
500
1000
1500 2000 2500 Frequency (cm-1)
3000
3500
(d) ref torsion
250 200 150 100 50 0
200
400
600
800
1000
time (fs)
FIG. 2. Analysis of two representative ADMP trajectories starting at the minimized B3LYP/LANL2DZ structure (ref) and from a structure displaced along the C2 -C1 -S-Ag dihedral bend (torsion, E from the minimum = 2.6 × 10−2 kcal mol−1 ) at t = 0. (a) Changes in the total and potential energy of the system on a 1 ps timescale. (b) Fourier transform of the nuclear velocity autocorrelation functions. (c) and (d): Time evolution of the C1 -S-Ag angle and the C2 -C1 -S-Ag dihedral. See text for more details.
E (displaced-minimized) = 2.6 × 10−2 kcal mol−1 at t = 0. In both cases, the nuclear velocities are randomly partitioned so that the total nuclear kinetic energy of the system is held at 0.1 Hartree. We illustrate the behavior of the respective Fourier transforms of the velocity autocorrelation functions (FT-VAC) derived from the selected AMDP trajectories in Figure 2(b), noting that the FT-VACs calculated from structures displaced with different increments along the C2 -C1 -SAg dihedral are nearly identical. We average the five FT-VACs obtained by translating the optimized structure by different increments along the dihedral bend, see Figures 1(b) and 2(b). Similarly, we compute an additional set of 5 trajectories, initiated by displacing the minimum along the aromatic C = C stretching vibration (ν 46 = 1611 cm−1 ). Detailed information about the calculations and their analysis can be found in the supplementary material.17 III. RESULTS AND DISCUSSION
The harmonic infrared and Raman spectra computed at the minimum energy conformation to the two FT-VACs are compared in Figure 1(b). The FT-VAC contains information on all the ro-vibrational density of states sampled throughout the trajectory, and only those modes which exhibit a change in dipole moment and polarizability would manifest in the infrared and Raman spectra, respectively. Distinct differences are observed between the FT-VAC computed from the reference trajectory and its analogue accelerated along the torsional mode, see Figure 2(b). Notable differences in the 1000–2000 cm−1 region of interest to our on-going SERS experiments are the relatively large amplitudes of the 1571 and 1078 cm−1 modes in the latter. The two vibrational states
are associated with the totally symmetric ring stretching and the C-S stretching modes, respectively. One way to rationalize the observed differences between the two FT-VACs is to compare the changes in key geometrical parameters along the computed ADMP trajectories. We analyze the C1 -S-Ag angle and C2 -C1 -S-Ag dihedral, with a premise that motion in the space spanned by these two parameters controls the coupling between molecule and finite metal cluster, see Figures 2(c) and 2(d).10 The angle oscillates by the same extent around its optimized value of ∼100◦ in the two trajectories. On the other hand, initial “acceleration” about the torsional motion causes an almost complete rotation about the C-S bond on a 1 ps timescale in the torsion-mediated trajectory, not observed in the reference trajectory which gently oscillates about its optimized value of ∼77◦ . It appears that the low frequency torsional motion which controls the vibrational density of states in our model system is strongly coupled to higher frequency vibrations in the 1000–2000 cm−1 region of interest. Trajectories initiated by displacing the molecule along the aromatic C = C stretching vibration further validate this premise.17 They reveal that displacing the optimized structure along the stretch also leads to a complete rotation about the C-S bond throughout the course of the 1 ps trajectories. How facile torsional motion would in turn be expected to affect the computed UV-Vis, IR, and Raman spectra is addressed next. Figure 3(a) qualitatively illustrates the change in the Highest Occupied Molecular Orbital (HOMO) electron density along the above-mentioned torsion-mediated ADMP trajectory. Notice the sloshing of the HOMO electron density in/out of the π -orbitals of the aromatic ring as a function of propagation time. Figure 3(b) illustrates the drastic change in the HOMO-LUMO gap as a function of the C2 -C1 -S-Ag
174303-4
El-Khoury, Bylaska, and Hess
J. Chem. Phys. 139, 174303 (2013)
2.00
TABLE II. A comparison of the B3LYP/LANL2DZ Raman active vibrations computed at the minimized structure and derived from ADMP trajectories (average of 5 ν9-displaced FT-VACs, see Figures 1 and 2) to their experimental (Ref. 14) and computed (Ref. 10) analogues. Values are given in cm−1 .
1.75
Source
2.50
HOMO-LUMO Gap (eV)
(a)
2.25
1.50
(b) 1.25 50 75 100 125 150 175 200 225 250 275 300 325 7
Diplole Moment (Debye)
6
ω2
ω3
ω4
ω5
ω6
1576 1550 1566 1611 1571
1076 1063 1070 1089 1078
1027 1018 1022 1038 1028
1003 997 1000 1002 ...
695 705 717 689 708
422 415 417 419 431
5 4 3 2 x
1
z
(c)
0 50 75 100 125 150 175 200 225 250 275 300 325 410 400
Polarizability (Bohr3)
Expt.–Ref. 14 Reference 10–I Reference 10–II This work–minimum This work–ADMP
ω1
390 380 370 50
isotropic xz
40 30
(d)
20
50 75 100 125 150 175 200 225 250 275 300 325 Dihedral Angle (Degrees)
FIG. 3. Analysis of the temporal evolution of the electronic structure ((a) and (b)), x- and z-components of the electric dipole moment (μx and μz , (c)), isotropic polarizability, and xz-component of the polarizability tensor (α xz , (d)) along the computed ADMP trajectory accelerated along the torsional motion, and as shown in Figures 1 and 2.
dihedral angle along torsion-mediated trajectory, reaching a minimum value of 1.4 eV at ∼167◦ , where the π -framework of the benzyl moiety is decoupled from the Ag cluster. Metalto-molecule charge transfer accounts for the observed destabilization/stabilization of the HOMO/LUMO in the latter, a detailed analysis of which goes beyond the realm of the present work, see the supplementary material.17 On the other hand, maximum gap values of 2.34 and 2.38 eV are attained at < C2 -C1 -S-Ag dihedral angle values of 77◦ (t = 200 fs, structure is shown in Figure 3(a)) and 254◦ , respectively, where the interaction between the π ring and the cluster is strongest. An angular difference of 177◦ and the shape of the trace in Figure 3(b) both reflect the symmetry in rotating about the C-S bond. Recognizing the former, we monitor the changes in the x and z components of the dipole (μx and μz , see Figure 3(c)) as well as the change in the isotropic and xz (α xz ) components of the molecular polarizability tensors (see Figure 3(d)) as a function of the dihedral angle (synonymous with time in this case). An inverse correlation between μx and μz is observed, reflecting a change in the dipole direction, which can also be visualized by inspecting the evolution of the blue arrows superimposed on the molecular structure and electron densities in Figure 3(a). Changes in the tensor property are less trivial to readily interpret and visualize in the ab-
sence of analyzing incident and scattered radiation fields (e.g., in Raman spectroscopy). That said, Figure 3(d) illustrates the inverse correlation between the α xz and isotropic components of the molecular polarizability tensor along the computed trajectory. This study suggests that sampling a single conformation does not yield an accurate description of the electronic structure and vibrational properties of PhS-Ag7 , and likely of other SERS model systems. This shortcoming can be circumvented by tackling the problem in the time domain as illustrated herein. In Table II, we demonstrate that although our infrared and Raman simulations which are computed at the minimized vertex-bound structure do not match well with their experimental counterparts as compared to the superior level of theory used in the previous work,10 the vibrational structure derived from the time-domain simulations in this work are well aligned with the reported experimental values.14 We acknowledge that the agreement between our simulations and previous experiments may be fortuitous, given the small basis set employed in the current report. Taken at face value, our observations underscore the important role played by large amplitude motions in controlling moleculemetal coupling in experiment, as captured by our choice of initial conditions in the reported AIMD simulations which also sample the anharmonic potential energy surface. Overall, the significant changes in the dipole, polarizability, and HOMO-LUMO energy gap as a function of the dihedral angle all suggest that the predicted properties will be significantly affected by the relative orientation of the ring with respect to the Ag cluster.
IV. CONCLUSIONS
In conclusion, our present study suggests that AIMD calculations provide invaluable insights into chemical bonding effects in surface-enhanced vibrational spectroscopy. In this regard, more systematics are warranted in follow-up reports. Specifically, longer trajectories would be needed to better account for the effects of low-frequency vibrational modes in the spectrum, as the number of configurations which contribute to the spectrum grows exponentially as a result of the low frequency/large amplitude modes present in model systems similar to that considered in this study. Herein, we identify the important role played by the torsional motion — which controls the coupling between the silver cluster and the
174303-5
El-Khoury, Bylaska, and Hess
aromatic ring of the chemisorbed molecule — in inducing significant changes to the predicted dipoles, polarizabilities, and orbital energies/densities. On-going works from our group aim at correlating spectral fluctuations observed in SERS experiments with molecular motion, in light of the conclusions drawn from this study.
ACKNOWLEDGMENTS
W.P.H. acknowledges support from the US Department of Energy (DOE), Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences and Biosciences. P.Z.E. acknowledges support from the Laboratory Directed Research and Development Program through a Linus Pauling Fellowship at Pacific Northwest National Laboratory (PNNL), an allocation of computing time from the National Science Foundation (TG-CHE130003), and the use of the Extreme Science and Engineering Discovery Environment. This work was performed using EMSL, a national scientific user facility sponsored by the DOE’s Office of Biological and Environmental Research and located at PNNL. PNNL is a multiprogram national laboratory operated for DOE by Battelle.
J. Chem. Phys. 139, 174303 (2013) 1 M.
Fleischmann, P. J. Hendra, and A. J. McQuillan, Chem. Phys. Lett. 26, 163 (1974). 2 D. L. Jeanmaire and R. P. Van Duyne, J. Electroanal. Chem. 84, 1 (1977). 3 M. Moskovits, J. Chem. Phys. 69, 4159 (1978). 4 A. Otto, Surf. Sci. 75, L392 (1978). 5 K. Kneipp, Y. Wang, H. Kneipp, L. Perelman, I. Itzkan, R. Dasari, and M. Feld, Phys. Rev. Lett. 78, 1667 (1997). 6 S. Nie and S. R. Emory, Science 275, 1102 (1997). 7 E. Le Ru and P. G. Etchegoin, Annu. Rev. Phys. Chem. 63, 65 (2012). 8 L. Jensen, C. M. Aikens, and G. C. Schatz, Chem. Soc. Rev. 37, 1061 (2008). 9 L. Zhao, L. Jensen, and G. C. Schatz, J. Am. Chem. Soc. 128, 2911 (2006). 10 S. K. Saikin, R. Olivared-Amaya, D. Rappoport, M. Stopa, and A. Aspuru-Gizik, Phys. Chem. Chem. Phys. 11, 9401 (2009). 11 H. B. Schlegel, J. M. Millam, S. S. Iyengar, G. A. Voth, G. E. Scuseria, A. D. Daniels, and M. J. Frisch, J. Chem. Phys. 114, 9758 (2001). 12 S. S. Iyengar, H. B. Schlegel, J. M. Millam, G. A. Voth, G. E. Scuseria, and M. J. Frisch, J. Chem. Phys. 115, 10291 (2001). 13 H. B. Schlegel, S. S. Iyengar, X. Li, J. M. Millam, G. A. Voth, G. E. Scuseria, and M. J. Frisch, J. Chem. Phys. 117, 8694 (2002). 14 K. B. Biggs, J. P. Camden, J. N. Anker, and R. P. Van Duyne, J. Phys. Chem. A 113, 4581 (2009). 15 M. J. Frisch, G. W. Trucks, H. B. Schlegel et al., GAUSSIAN 09, C.01, Gaussian, Inc., Wallingford, CT, 2009. 16 G. H. Golub and C. F. van Loan, Matrix Computations (The Johns Hopkins University Press, Baltimore, 1996). 17 See supplementary material at http://dx.doi.org/10.1063/1.4827455 for a detailed analysis of the individual ADMP trajectories discussed in the main text and for an analysis of the HOMO and LUMO orbitals.
The Journal of Chemical Physics is copyrighted by the American Institute of Physics (AIP). Redistribution of journal material is subject to the AIP online journal license and/or AIP copyright. For more information, see http://ojps.aip.org/jcpo/jcpcr/jsp
Time domain simulations of chemical bonding effects in surface-enhanced spectroscopy Patrick Z. El-Khoury,a) Eric J. Bylaska, and Wayne P. Hessb) Physical Sciences Division, Pacific Northwest National Laboratory, P.O. Box 999, Richland, Washington 99352, USA
(Received 12 July 2013; accepted 16 October 2013; published online 4 November 2013) The atom-centered density-matrix propagation method is used to illustrate how time-dependent conformational changes affect the electronic structure and derived spectroscopic properties of a prototypical finite metal cluster-bound π -conjugated organic complex, Ag7 -benzenethiol. We establish that there is considerable conformational flexibility to the model structure, even at relatively low temperatures, which influences the predicted spectroscopic properties. Namely, the computed electron densities, dipoles, and polarizabilities are all dictated by torsional motion which controls the coupling between the π -framework of the chemisorbed molecular system and the cluster. © 2013 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4827455] I. INTRODUCTION
To account for the underlying dynamics in surfaceenhanced spectroscopy, it is essential to consider the distinct local environments in which different molecules reside. Electric field localization at defect sites in metal substrates induces a selective enhancement of the optical response of nearby molecules. This is a direct consequence of the inherent heterogeneity of typical metallic terrains over which molecular reporters are physi- or chemi-sorbed. In this regard, early works in the area of surface-enhanced spectroscopy1–4 inspired several ensuing studies,5–7 in which engineered hotspots were found to afford single molecule detection sensitivity. Understanding molecule-metal coupling in hybrid molecularmetallic constructs is a real challenge, both from the theoretical and experimental standpoints. It is a problem which lies at the interface of classical and quantum theories — each of which is individually tailored to solely describe either molecule or metal. A practical approach to rationalize chemical bonding effects in SERS is to simulate the local environment of a molecule using static tools of density functional theory (DFT) and its time-dependent variant (TD-DFT).8 In a pioneering work by Schatz and co-workers, the localization of electromagnetic fields at hotspots within a metallic substrate is identified as one of three cooperative mechanisms in surface-enhanced Raman spectroscopy (SERS).9 The other two are chemical in nature, namely (i) changes in the static polarizability of a molecule upon adsorption, and (ii) resonant enhancement due to metal-molecule charge transfer transitions. Of direct relevance to our present report is another early study by Aspuru-Guzik and co-workers, where the effects of chemical bonding on Raman scattering from benzenethiol (BT) chemisorbed on silver clusters was investigated using DFT calculations.10 Using different finite silver
a) [email protected] b) [email protected]
0021-9606/2013/139(17)/174303/5/$30.00
cluster models (Agn , n = 6–11), their work illustrated the need to consider various conformations of BT with respect to the finite cluster to correlate the metal-induced changes in the electronic structure of a molecule with SERS enhancement. Of the different effects governing normal (off-resonance) SERS enhancement, three were highlighted, namely (i) the relative orientation of the π framework with respect to the cluster, (ii) the local symmetry of the benzene ring, and (iii) the proximity of the different vibrational modes with respect to the binding sites. Our current work is driven by these observations and general interest in SERS. Rather than considering different clusters to emulate the aforementioned effects using static tools of quantum chemical simulations, we perform ab initio molecular dynamics (AIMD) simulations at the DFT level to systematically explore operative mechanisms in hotspot spectroscopy using one model system. Beyond the static picture, i.e., computing the optical properties starting at energy-minimized structures, the effect of time-dependent conformational and structural rearrangement in finite metal cluster-bound systems on their spectroscopic properties remains uninvestigated. Recognizing that the latter can at the very least be driven by thermal population under ambient laboratory conditions, we undertake the challenge of correlating the temporal evolution of molecular structure on a picosecond timescale with the derived electronic and vibrational structures in a prototypical model system, namely Ag7 -bound BT (PhS-Ag7 ). To this end, we perform AIMD simulations using the atom-centered density-matrix propagation (ADMP) method of Schlegel et al.11–13 The PhS-Ag7 system was chosen as a prototypical system because it has been well-characterized in previous theoretical10 and experimental works.14 This work extends previous theoretical studies by correlating time-dependent structural evolution with spectroscopic properties. We establish that conformational flexibility to our model structure dictates the computed electron densities, orbital energies (HOMO-LUMO gaps), dipole moments, and polarizabilities.
139, 174303-1
© 2013 AIP Publishing LLC
174303-2
El-Khoury, Bylaska, and Hess
J. Chem. Phys. 139, 174303 (2013)
II. COMPUTATIONAL METHODS
The B3LYP exchange-correlation functional and the LANL2DZ basis set are used throughout this study. Of the two possible isomers of PhS-Ag7 , we select the higher energy conformer (the vertex bound PhS-Ag7 structure shown in Figure 1) as a starting point for geometry optimization and frequency analysis. In principle, our calculations would allow for the interconversion between the two isomer forms where enough initial kinetic energy is partitioned into the isomerization coordinate. The harmonic IR and Raman spectra are represented as sums of Lorentzian functions that were empirically broadened by 15 cm−1 . In this concept report, the reduced basis set (266 basis functions) was chosen to decrease the computational burden of our rather resourceexhaustive simulations. However, in comparing our computed B3LYP/LANL2DZ properties at the minimum to their counterparts at the PBE0/def2-TZVP level of theory,10 it appears that the reduced representation is reasonable for the purpose of our current report, see Table I. Further validation of the level of theory used in the ADMP scheme follows. In Table I, we estimate the binding energy (Eb ) between BT and the Ag7 cluster is according to Eb (P hS − Ag7 ) = E(Ag7 ) + E(P hSH ) − E(P hS − Ag7 ) − 0.5E(H2 ).
(1)
We neglect basis set superposition errors and find that the computed Eb is reasonable, falling within the range of binding energies previously reported for BT bound to different sizes of silver clusters.10 The GAUSSIAN 0915 implementation of the ADMP methodology was used to carry out the AIMD simulations. The method is well-described in the original publications.11–13 Briefly, these are Car-Parrinello type AIMD simulations, carried out by propagating the density matrix and performed using atom centered Gaussians. The density matrices were obtained at each time step (t = 0.2 fs) from fully converged
TABLE I. A comparison of (i) selected structural parameters (bond lengths, angles in Ångstroms, degrees), (ii) binding energies (kcal mol−1 ), (iii) isotropic polarizabilities (a.u.), and (iv) HOMO-LUMO gaps (eV) in the PhS-Ag7 model computed in this work with their counterparts from Ref. 10. C2 -Ag denoted the shortest distance between the aromatic ring (a carbon atom adjacent to the C1 -S bond) and the cluster. The binding energy is calculated according to Eq. (1) in the main text. Source
Ag-S
C-S
This work 2.528 1.844 Reference 10–I 2.421 1.757 Reference 10–II 2.521/ 1.768 2.589
< C-S-Ag C2 -Ag 99.795 96.606 103.673
2.719 2.669 2.629
Eb
α iso
Egap
22.05 371.5 2.34 17.71 365.2 2.80 20.67 374.9 2.54
SCF calculations. The relatively small time step is chosen such that the total energy was converged to within ∼10−4 hartree. In these calculations, we use the same level of theory used for geometry optimization and harmonic frequency calculations, and employ the Cholesky decomposition scheme.16 We began by computing ADMP trajectories starting from the minimized structure. The initial nuclear velocities were randomly partitioned so that the total nuclear kinetic energies were 0.01, 0.05, 0.1, 0.15, and 0.2 hartree for five different trajectories of over 1 ps in total time duration. Next, to study the effects of large amplitude motion on the derived electronic and vibrational structures in the system, we identify a normal mode (ν 9 = 70 cm−1 ) of the vertex-bound structure that is best described as torsional motion along the C2 -C1 -S-Ag dihedral, see Figure 1(a). We translate the optimized structure by different increments along ν 9 , and use the resulting geometries as starting structures for another set of ADMP trajectories. Of the ADMP simulations performed, we analyze two representative trajectories in Figures 1 and 2. The first is initiated from the vertex-bound structure (henceforth designated as the reference trajectory). The starting structure for the second trajectory was obtained by displacing the same structure along the C2 -C1 -S-Ag dihedral bend, such that
0.6
(b)
Normalized Intensities (a.u.)
Relative energy (kcal mol-1)
(a) 0.5
0.4
0.3
0.2
0.1
0.0 0.0
0.1
0.2
0.3
Displacement along 9 (arbitary units)
0.4
0
500
1000 1500 2000 2500 3000 3500 Frequency (cm-1)
FIG. 1. (a) Change in energy for translating the minimized structure (0,0) by different increments along the C2 -C1 -S-Ag dihedral bend (torsion), a vector representation of which is shown in the inset. (b) Comparison of the Infrared and Raman spectra computed at the minimum to the Fourier transform of the nuclear velocity autocorrelation functions derived from ADMP simulations. The black trace (bottom) corresponds to the vibrational spectrum derived from the reference ADMP trajectory. Five frequency traces obtained by translating the minimum energy structure by different increments along the C2 -C1 -S-Ag dihedral bend at t = 0 were averaged to obtain the spectrum shown (dark yellow trace, second from the bottom).
J. Chem. Phys. 139, 174303 (2013)
Energy (Hartree)
-1262.26 TE - ref
TE - torsion
-1262.28 -1262.30 -1262.32 -1262.34
PE - ref
PE - torsion
(a)
-1262.36 0
200
400 600 time (fs)
800
C-S-Ag Angle (Degrees)
120 110 100 90 80
(c)
70 0
200
400
600
800
1.0
1000
time (fs)
(b)
0.8
FFT - ref FFT - torsion
0.6 0.4 0.2 0.0
1000
130
ref torsion
Normalized Amplitude (a.u.)
El-Khoury, Bylaska, and Hess
C-C-S-Ag Dihedral (Degrees)
174303-3
0
300
500
1000
1500 2000 2500 Frequency (cm-1)
3000
3500
(d) ref torsion
250 200 150 100 50 0
200
400
600
800
1000
time (fs)
FIG. 2. Analysis of two representative ADMP trajectories starting at the minimized B3LYP/LANL2DZ structure (ref) and from a structure displaced along the C2 -C1 -S-Ag dihedral bend (torsion, E from the minimum = 2.6 × 10−2 kcal mol−1 ) at t = 0. (a) Changes in the total and potential energy of the system on a 1 ps timescale. (b) Fourier transform of the nuclear velocity autocorrelation functions. (c) and (d): Time evolution of the C1 -S-Ag angle and the C2 -C1 -S-Ag dihedral. See text for more details.
E (displaced-minimized) = 2.6 × 10−2 kcal mol−1 at t = 0. In both cases, the nuclear velocities are randomly partitioned so that the total nuclear kinetic energy of the system is held at 0.1 Hartree. We illustrate the behavior of the respective Fourier transforms of the velocity autocorrelation functions (FT-VAC) derived from the selected AMDP trajectories in Figure 2(b), noting that the FT-VACs calculated from structures displaced with different increments along the C2 -C1 -SAg dihedral are nearly identical. We average the five FT-VACs obtained by translating the optimized structure by different increments along the dihedral bend, see Figures 1(b) and 2(b). Similarly, we compute an additional set of 5 trajectories, initiated by displacing the minimum along the aromatic C = C stretching vibration (ν 46 = 1611 cm−1 ). Detailed information about the calculations and their analysis can be found in the supplementary material.17 III. RESULTS AND DISCUSSION
The harmonic infrared and Raman spectra computed at the minimum energy conformation to the two FT-VACs are compared in Figure 1(b). The FT-VAC contains information on all the ro-vibrational density of states sampled throughout the trajectory, and only those modes which exhibit a change in dipole moment and polarizability would manifest in the infrared and Raman spectra, respectively. Distinct differences are observed between the FT-VAC computed from the reference trajectory and its analogue accelerated along the torsional mode, see Figure 2(b). Notable differences in the 1000–2000 cm−1 region of interest to our on-going SERS experiments are the relatively large amplitudes of the 1571 and 1078 cm−1 modes in the latter. The two vibrational states
are associated with the totally symmetric ring stretching and the C-S stretching modes, respectively. One way to rationalize the observed differences between the two FT-VACs is to compare the changes in key geometrical parameters along the computed ADMP trajectories. We analyze the C1 -S-Ag angle and C2 -C1 -S-Ag dihedral, with a premise that motion in the space spanned by these two parameters controls the coupling between molecule and finite metal cluster, see Figures 2(c) and 2(d).10 The angle oscillates by the same extent around its optimized value of ∼100◦ in the two trajectories. On the other hand, initial “acceleration” about the torsional motion causes an almost complete rotation about the C-S bond on a 1 ps timescale in the torsion-mediated trajectory, not observed in the reference trajectory which gently oscillates about its optimized value of ∼77◦ . It appears that the low frequency torsional motion which controls the vibrational density of states in our model system is strongly coupled to higher frequency vibrations in the 1000–2000 cm−1 region of interest. Trajectories initiated by displacing the molecule along the aromatic C = C stretching vibration further validate this premise.17 They reveal that displacing the optimized structure along the stretch also leads to a complete rotation about the C-S bond throughout the course of the 1 ps trajectories. How facile torsional motion would in turn be expected to affect the computed UV-Vis, IR, and Raman spectra is addressed next. Figure 3(a) qualitatively illustrates the change in the Highest Occupied Molecular Orbital (HOMO) electron density along the above-mentioned torsion-mediated ADMP trajectory. Notice the sloshing of the HOMO electron density in/out of the π -orbitals of the aromatic ring as a function of propagation time. Figure 3(b) illustrates the drastic change in the HOMO-LUMO gap as a function of the C2 -C1 -S-Ag
174303-4
El-Khoury, Bylaska, and Hess
J. Chem. Phys. 139, 174303 (2013)
2.00
TABLE II. A comparison of the B3LYP/LANL2DZ Raman active vibrations computed at the minimized structure and derived from ADMP trajectories (average of 5 ν9-displaced FT-VACs, see Figures 1 and 2) to their experimental (Ref. 14) and computed (Ref. 10) analogues. Values are given in cm−1 .
1.75
Source
2.50
HOMO-LUMO Gap (eV)
(a)
2.25
1.50
(b) 1.25 50 75 100 125 150 175 200 225 250 275 300 325 7
Diplole Moment (Debye)
6
ω2
ω3
ω4
ω5
ω6
1576 1550 1566 1611 1571
1076 1063 1070 1089 1078
1027 1018 1022 1038 1028
1003 997 1000 1002 ...
695 705 717 689 708
422 415 417 419 431
5 4 3 2 x
1
z
(c)
0 50 75 100 125 150 175 200 225 250 275 300 325 410 400
Polarizability (Bohr3)
Expt.–Ref. 14 Reference 10–I Reference 10–II This work–minimum This work–ADMP
ω1
390 380 370 50
isotropic xz
40 30
(d)
20
50 75 100 125 150 175 200 225 250 275 300 325 Dihedral Angle (Degrees)
FIG. 3. Analysis of the temporal evolution of the electronic structure ((a) and (b)), x- and z-components of the electric dipole moment (μx and μz , (c)), isotropic polarizability, and xz-component of the polarizability tensor (α xz , (d)) along the computed ADMP trajectory accelerated along the torsional motion, and as shown in Figures 1 and 2.
dihedral angle along torsion-mediated trajectory, reaching a minimum value of 1.4 eV at ∼167◦ , where the π -framework of the benzyl moiety is decoupled from the Ag cluster. Metalto-molecule charge transfer accounts for the observed destabilization/stabilization of the HOMO/LUMO in the latter, a detailed analysis of which goes beyond the realm of the present work, see the supplementary material.17 On the other hand, maximum gap values of 2.34 and 2.38 eV are attained at < C2 -C1 -S-Ag dihedral angle values of 77◦ (t = 200 fs, structure is shown in Figure 3(a)) and 254◦ , respectively, where the interaction between the π ring and the cluster is strongest. An angular difference of 177◦ and the shape of the trace in Figure 3(b) both reflect the symmetry in rotating about the C-S bond. Recognizing the former, we monitor the changes in the x and z components of the dipole (μx and μz , see Figure 3(c)) as well as the change in the isotropic and xz (α xz ) components of the molecular polarizability tensors (see Figure 3(d)) as a function of the dihedral angle (synonymous with time in this case). An inverse correlation between μx and μz is observed, reflecting a change in the dipole direction, which can also be visualized by inspecting the evolution of the blue arrows superimposed on the molecular structure and electron densities in Figure 3(a). Changes in the tensor property are less trivial to readily interpret and visualize in the ab-
sence of analyzing incident and scattered radiation fields (e.g., in Raman spectroscopy). That said, Figure 3(d) illustrates the inverse correlation between the α xz and isotropic components of the molecular polarizability tensor along the computed trajectory. This study suggests that sampling a single conformation does not yield an accurate description of the electronic structure and vibrational properties of PhS-Ag7 , and likely of other SERS model systems. This shortcoming can be circumvented by tackling the problem in the time domain as illustrated herein. In Table II, we demonstrate that although our infrared and Raman simulations which are computed at the minimized vertex-bound structure do not match well with their experimental counterparts as compared to the superior level of theory used in the previous work,10 the vibrational structure derived from the time-domain simulations in this work are well aligned with the reported experimental values.14 We acknowledge that the agreement between our simulations and previous experiments may be fortuitous, given the small basis set employed in the current report. Taken at face value, our observations underscore the important role played by large amplitude motions in controlling moleculemetal coupling in experiment, as captured by our choice of initial conditions in the reported AIMD simulations which also sample the anharmonic potential energy surface. Overall, the significant changes in the dipole, polarizability, and HOMO-LUMO energy gap as a function of the dihedral angle all suggest that the predicted properties will be significantly affected by the relative orientation of the ring with respect to the Ag cluster.
IV. CONCLUSIONS
In conclusion, our present study suggests that AIMD calculations provide invaluable insights into chemical bonding effects in surface-enhanced vibrational spectroscopy. In this regard, more systematics are warranted in follow-up reports. Specifically, longer trajectories would be needed to better account for the effects of low-frequency vibrational modes in the spectrum, as the number of configurations which contribute to the spectrum grows exponentially as a result of the low frequency/large amplitude modes present in model systems similar to that considered in this study. Herein, we identify the important role played by the torsional motion — which controls the coupling between the silver cluster and the
174303-5
El-Khoury, Bylaska, and Hess
aromatic ring of the chemisorbed molecule — in inducing significant changes to the predicted dipoles, polarizabilities, and orbital energies/densities. On-going works from our group aim at correlating spectral fluctuations observed in SERS experiments with molecular motion, in light of the conclusions drawn from this study.
ACKNOWLEDGMENTS
W.P.H. acknowledges support from the US Department of Energy (DOE), Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences and Biosciences. P.Z.E. acknowledges support from the Laboratory Directed Research and Development Program through a Linus Pauling Fellowship at Pacific Northwest National Laboratory (PNNL), an allocation of computing time from the National Science Foundation (TG-CHE130003), and the use of the Extreme Science and Engineering Discovery Environment. This work was performed using EMSL, a national scientific user facility sponsored by the DOE’s Office of Biological and Environmental Research and located at PNNL. PNNL is a multiprogram national laboratory operated for DOE by Battelle.
J. Chem. Phys. 139, 174303 (2013) 1 M.
Fleischmann, P. J. Hendra, and A. J. McQuillan, Chem. Phys. Lett. 26, 163 (1974). 2 D. L. Jeanmaire and R. P. Van Duyne, J. Electroanal. Chem. 84, 1 (1977). 3 M. Moskovits, J. Chem. Phys. 69, 4159 (1978). 4 A. Otto, Surf. Sci. 75, L392 (1978). 5 K. Kneipp, Y. Wang, H. Kneipp, L. Perelman, I. Itzkan, R. Dasari, and M. Feld, Phys. Rev. Lett. 78, 1667 (1997). 6 S. Nie and S. R. Emory, Science 275, 1102 (1997). 7 E. Le Ru and P. G. Etchegoin, Annu. Rev. Phys. Chem. 63, 65 (2012). 8 L. Jensen, C. M. Aikens, and G. C. Schatz, Chem. Soc. Rev. 37, 1061 (2008). 9 L. Zhao, L. Jensen, and G. C. Schatz, J. Am. Chem. Soc. 128, 2911 (2006). 10 S. K. Saikin, R. Olivared-Amaya, D. Rappoport, M. Stopa, and A. Aspuru-Gizik, Phys. Chem. Chem. Phys. 11, 9401 (2009). 11 H. B. Schlegel, J. M. Millam, S. S. Iyengar, G. A. Voth, G. E. Scuseria, A. D. Daniels, and M. J. Frisch, J. Chem. Phys. 114, 9758 (2001). 12 S. S. Iyengar, H. B. Schlegel, J. M. Millam, G. A. Voth, G. E. Scuseria, and M. J. Frisch, J. Chem. Phys. 115, 10291 (2001). 13 H. B. Schlegel, S. S. Iyengar, X. Li, J. M. Millam, G. A. Voth, G. E. Scuseria, and M. J. Frisch, J. Chem. Phys. 117, 8694 (2002). 14 K. B. Biggs, J. P. Camden, J. N. Anker, and R. P. Van Duyne, J. Phys. Chem. A 113, 4581 (2009). 15 M. J. Frisch, G. W. Trucks, H. B. Schlegel et al., GAUSSIAN 09, C.01, Gaussian, Inc., Wallingford, CT, 2009. 16 G. H. Golub and C. F. van Loan, Matrix Computations (The Johns Hopkins University Press, Baltimore, 1996). 17 See supplementary material at http://dx.doi.org/10.1063/1.4827455 for a detailed analysis of the individual ADMP trajectories discussed in the main text and for an analysis of the HOMO and LUMO orbitals.
The Journal of Chemical Physics is copyrighted by the American Institute of Physics (AIP). Redistribution of journal material is subject to the AIP online journal license and/or AIP copyright. For more information, see http://ojps.aip.org/jcpo/jcpcr/jsp
