Theory for Modeling of High Resolution Resonant and Nonresonant

Sep 20, 2016 - Journal of Chemical Education · Journal of Chemical Information .... Department of Theoretical Chemistry and Biology, School of Biotech...
1 downloads 0 Views 2MB Size
Article pubs.acs.org/JCTC

Theory for Modeling of High Resolution Resonant and Nonresonant Raman Images Sai Duan,†,∥ Guangjun Tian,‡,∥ and Yi Luo*,§,† †

Department of Theoretical Chemistry and Biology, School of Biotechnology, Royal Institute of Technology, S-106 91 Stockholm, Sweden ‡ College of Science, Yanshan University, Qinhuangdao 066004, China § Hefei National Laboratory for Physical Science at the Microscale, Department of Chemical Physics, School of Chemistry and Materials Science, University of Science and Technology of China, Hefei, 230026 Anhui, China ABSTRACT: Tip-enhanced Raman imaging is capable of resolving the inner structure of a single molecule owing to the generation of highly localized nanocavity plasmon. Here we present a general theory and detailed computational methodology to fully describe resonant and nonresonant Raman scattering under the localized plasmonic field. We use an allylcarbinol molecule adsorbed on the gold surface as a model system to illustrate different effects on the Raman images. It is found that the ability of distinguishing an individual vibration mode is highly limited under the resonant condition due to the dominant contribution from the Franck−Condon term and the mode-independent component of the Herzberg−Teller term. The nonresonant Raman images of the single molecule are vibrationally distinguishable and present the vibrational motion of the corresponding vibrational modes in real space. Furthermore, the calculated results confirm that nonlinear optical effects can further improve the resolution of the images. The theoretical and computational methods presented here provide the basic tools to model high resolution Raman images at the single molecular level.

1. INTRODUCTION Optical imaging techniques are widely used in surface science, molecular nanotechnology, and biotechnology, which are motivated by the spirit of “seeing is believing”. However, because of the diffraction limit, the resolution of conventional far field optical images is quite low,1 even with the help of the super-resolution fluorescence techniques.2−4 In this context, by precisely controlling the localized surface plasmon that has huge near-field enhancement factor for incident electromagnetic (EM) field, tip-enhanced Raman spectroscopy (TERS)5−8 has become a promising practical tool to break down the diffraction limit for the optical imaging of adsorbates.9−16 With the developments of the TERS technique, the resolution of the Raman images has been greatly improved over the past few years.9−16 For instance, in a recent experiment, the resolution of the Raman images has amazingly archived the sub-nanometer level under ultrahigh vacuum and low temperature conditions,15 which opens the door for the optical imaging of inner molecular structures. It is anticipated that, in the superhigh resolution Raman images, the spatially confined plasmon (SCP) between the TERS tip and substrate plays the decisive role.17 Due to the fact that the size of SCP is comparable to the molecular size, the widely used uniform treatment for the interactions between EM field and adsorbates as in conventional response theory18−25 becomes inappropriate. 26 In this case, the interaction Hamiltonian is dependent on the location as well as the © 2016 American Chemical Society

distribution of the SCP. In order to theoretically simulate the Raman images, we have proposed a theoretical framework to take into account the spatial distribution of the SCP in the calculation of the induced dipole moment.26 It was found that when the full width at half-maximum (fwhm) of the SCP was set to 20 Å in the scanning plane, the theoretical result can reproduce the experimental observations very well, highlighting the importance of the locality of the SCP. This predicted size of the SCP is consistent with the advanced EM simulations27−29 and further confirmed by recent experiments.30 It is noted that our method could be integrated with the electrodynamics simulations by the coordination-dependent discrete interaction model31 as well as with the multiscale Maxwell−Schrödinger approaches.32−34 We also found that the similarity of the Raman images from different vibrational frequencies15 is due to the resonant conditions used in the experiments.26 Under the resonant conditions with the strongly dipole-allowed transition, the Raman intensity is dominated by the Franck−Condon (FC) term and vibrational integrals are independent from the position of the tip, which leads to vibrationally insensitive images.35 This result may hinder the applications of the Raman images in distinguishing vibrational modes. In principle, the contribution of the Herzberg−Teller (HT) term should be Received: June 9, 2016 Published: September 20, 2016 4986

DOI: 10.1021/acs.jctc.6b00592 J. Chem. Theory Comput. 2016, 12, 4986−4995

Article

Journal of Chemical Theory and Computation emerged to obtain the vibrationally dependent Raman images.35,36 In this context, the nonresonant Raman processes would be a natural choice. Our calculation results of water clusters adsorbed on the Au(111) surface demonstrated that nonresonant Raman images are indeed vibrationally resolved and moreover resemble the atomic movement of corresponding vibrational modes.36 Thus, the nonresonant Raman images can be used as a practical tool to visualize molecular vibrational modes in real space.36 We should emphasize that for weakly dipole-allowed transitions, the contribution of HT term would also be promoted.37,38 In the present work, we comprehensively investigate resonant and nonresonant Raman images of allylcarbinol (AC) which has three different functional groups, namely, one vinyl group that has a CC double bound, two methanediyl groups (CH2), and one O−H bond. Our calculation results show that the resonant Raman images are similar for some individual normal modes, while, the nonresonant images resemble the corresponding vibrational modes very well. The nonlinear optical effects are also discussed for the case of nonresonant Raman images.

ℏω Mi g(r)aê − ıωt + H.c. 2ϵ0V

Γ iabs →f =



μ0 ∂(ωμ) |τ |2 ⎤ 2 2 ⎥Mi |g| = V ϵ0 ∂ω |μ|2 ⎦

where ϵ is the unit vector of electric polarization, g(r) is the corresponding scale function of g(r), and we have the definition of

(7)

Here n is the photon number. By definition, we find that is the electric field amplitude of the incident light. We should stress that the current treatment is in the same framework of quantization for the propagated surface plasmons as proposed by Archambault et al.42 and the interaction Hamiltonian is consistent with the strength matrix element used in the theory of femtosecond stimulated Raman spectroscopy proposed by Lee and co-workers.21 Compared with the conventional response theory, to take the effects of spatial confinement for the SCP into account, the amplitude distribution function g should be kept inside in the integrals, which is the origination of Raman images with superhigh spatial resolution. It is nice to remark that, in the condition of Mi = 1 and g(r) = 1, eq 6 would naturally return to the expression in conventional response theory for uniform EM field. It is well-known that the linear Raman process consists of the absorption and the spontaneous emission, which can be generally described by Albrecht’s theory.18,24 For Raman images, because the emission signals can be detected in the far field, the matrix element related to the spontaneous emission is unchanged and only the matrix element related to absorption is modified by the distribution of the plasmonic field. For restricting resonant conditions, a specific excited state is usually selected. Therefore, in analogy to Albrecht’s theory, the polarizability for linear Raman scattering that takes the spatial distribution of SCP into account could be calculated straightforwardly as26 E0i

(1)

αpq = A pq + Bpq

(2)

(8)

where A pq =

∞ FP ⎡ ⟨v f |v r ⟩⟨v r |v i⟩ ⎤ ⎢⟨ψ |p ̂|ψ ⟩⟨ψ |qĝ (r)|ψ ⟩ ∑ ⎥ + NRT g r r g ℏ ⎢⎣ ωervr : e g vi − ω − ıΓ ⎥⎦ vr = 0

(9) Bpq =

(3)

It should be emphasized that here μ̂ contains all one-body electron position operators in the system, i.e., μ̂ = e∑ir̂i, where e is the charge of the electron. Substituting eq 1 into eq 3, the interaction Hamiltonian can be rewritten as a harmonic perturbation /̂ ′(r, t ) = v(̂ r)e− ıωt + H.c.

nℏω 2ϵ0V

Ei0 =

where ϵ is the relative permittivity, μ0 (μ) is the permeability of free space (relative permeability), and τ is the ratio of magnetic induction to electric field. Furthermore, the nonrelativistic interaction Hamiltonian between electrons and electric field of SCP can still be written as43−45 /̂ ′ = −μ ̂·E(̂ r, t )

2π e 2|MiEi0|2 2 |ϵ ·⟨Ψf |rĝ (r)|Ψ⟩| δ(Ef − Ei + ℏω) i ℏ (6)

where ℏ is the reduced Planck constant, ω is the frequency of the incident laser, ϵ0 is the permittivity of free space, V is the system volume, â is the annihilation operator for photon, Mi is the corresponding enhancement factor, g is a vector function corresponding to the amplitude distribution of the electric field for the SCP, and H.c. is the hermitian conjugate. Notice that, for specific nanostructures, g can be determined by the classical Maxwell’s equations39 if the quantum correction is properly addressed under the tunneling conditions.40 According to the energy conservation law in classical electromagnetic theory for metallic nanostructures without considering losses,41,42 the following normalization condition should be fulfilled

∫ dr 12 ⎢⎣ ∂(∂ωωϵ) +

(5)

is time independent. Therefore, based on Fermi’s golden rule, the transition rate for one-photon absorption between two molecular states can be expressed as44

2. METHODOLOGY In principle, the SCP can be regarded as a highly confined EM field. In this context, the electric field operator that takes the spatial confinement of the SCP into account may be written as E(̂ r, t ) = ı

ℏω Mi r·g(r)a ̂ 2ϵ0V

v(̂ r) = −eı

FP ℏ

⎡ ∂⟨ψ |p ̂|ψ ⟩

∑⎢ k

⎢⎣

+⟨ψg |p ̂|ψr ⟩

g

r

∂Q k



⟨ψr|qĝ (r)|ψg ⟩

∂⟨ψr|qĝ (r)|ψg ⟩ ∂Q k



∑ vr = 0

∑ vr = 0

⟨v f |Q k|v r ⟩⟨v r |v i⟩ ωervr : e g vi − ω − ıΓ

⟨v |v ⟩⟨v r |Q k|v i⟩ ⎤ ⎥ + NRT ωervr : e g vi − ω − ıΓ ⎥⎦ f

r

(10)

Here A and B represent FC and HT terms, respectively, p and q are Cartesian coordinates, FP is the Purcell factor46 that accounts for the enhancement of spontaneous emission, |ψg⟩ and |ψr⟩ are the electronic ground and resonant excited states, Qk is the corresponding normal mode, |vi⟩ and |vf⟩ are the initial

(4)

where 4987

DOI: 10.1021/acs.jctc.6b00592 J. Chem. Theory Comput. 2016, 12, 4986−4995

Article

Journal of Chemical Theory and Computation and final vibrational states in |ψg⟩, |vr⟩ is the vibrational state in |ψr⟩, ωervr:egvi is the frequency difference between |ψr⟩|vr⟩ and |ψg⟩|vi⟩, Γ is the damping factor, and NRT is the nonresonant term. We should emphasize the multimode nature of the vibrational states. For instance, we have the relationship |vr⟩ = | vr1v2r ...vrN⟩, where N is the number of vibrational modes. Interestingly, because of the difference between the absorption and emission operator, here the conventional relationship between the one-photon absorption cross-section and the imaginary part of the diagonal matrix elements for dynamic polarizability47,48 is not fulfilled. One can always use the plasmon modified emission matrix element to resume the relationship. However, the corresponding process would require that the emission is a localized plasmon and cannot be observed in the far field. Therefore, we do not consider this process in detail in the current work. It is worth noting here that both the A term and the first term of the B term are vibrationally independent and these two terms are dominant respectively for strongly dipole-allowed and very weakly dipoleallowed transitions. Using eq 8, the amplitude of the induced dipole could be calculated by18,24,49

μ0 = Miα Ei0

Under nonresonant conditions, the ω + ıΓ term could be safely neglected (the zero-frequency limit). Furthermore, using the state-to-state mapping between Albrecht’s theory and the perturbation theory53 and summation of all excited states, the amplitude of the induced dipole can be calculated as36 μ0, k =

αeff, pq = 2 ∑ r

(12)

where c is the speed of light, ν̃s is the wavenumber of scattering, and Md is the directional radiation pattern factor owing to the corresponding nanostructures.25 It is noted that, in practical calculations, the Purcell factor could be calculated by the radiated power ratio between a point dipole with and without metallic nanostructures.50 Previous simulations show that, usually the Purcell factor is approximately equal to Mi2.51 As a result, the scattering intensity in eq 12 has the scaling of Mi4, which naturally resumes the well-known |E|4 scaling in SERS.52 Under nonresonant conditions, a sum over all excited states for the resonant expression is required. We should emphasize that the NRT could contribute in this case. Specifically, the |er⟩ component of polarizability also has two terms as (13)

where ∞ FP ⎡ ⟨v f |v r ⟩⟨v r |v i⟩ ⎢⟨ψ |p ̂|ψ ⟩⟨ψ |qĝ (r)|ψ ⟩ ∑ g r r g r r g i ℏ ⎢⎣ v r = 0 ω e v : e v − (ω + ı Γ)

r = A pq



+⟨ψg |qĝ (r)|ψr ⟩⟨ψr|p ̂|ψg ⟩

∑ r

v =0

FP

r Bpq =



⎡ ∂⟨ψ |p ̂|ψ ⟩

∑⎢ k

g

⎢⎣

+ ⟨ψg |p ̂|ψr ⟩ +

⎤ ⟨v f |v r ⟩⟨v r |v i⟩ ⎥ ωervr : e g v f + (ω + ıΓ) ⎥⎦

r

∂Q k



⟨ψr|qĝ (r)|ψg ⟩

∂⟨ψr|qĝ (r)|ψg ⟩ ∂Q k

∂⟨ψg |qĝ (r)|ψr ⟩ ∂Q k

+ ⟨ψg |qĝ (r)|ψr ⟩

⟨ψr|p ̂|ψg ⟩

∂⟨ψr|p ̂|ψg ⟩ ∂Q k



∑ vr = 0 ∞

∑ vr = 0 ∞

∑ r

v =0

ωervr : e g vi − (ω + ıΓ)

vr = 0 f

(14)

⟨v f |Q k|v r ⟩⟨v r |v i⟩

∑ r

⟨v |v ⟩⟨v r |Q k|v i⟩ ωervr : e g vi − (ω + ıΓ) ⟨v f |Q k|v r ⟩⟨v r |v i⟩ ωervr : e g v f + (ω + ıΓ) ⎤ ⎥ ωervr : e g v f + (ω + ıΓ) ⎥⎦

⟨ψg |p ̂|ψr ⟩⟨ψr|qĝ (r)|ψg ⟩ ΔErg

(17)

3. COMPUTATIONAL DETAILS To support the adsorbate, we used the Au(111) surface as the substrate. The Au(111) substrate is inertness for the adsorbates54 and could provide strong plasmonic coherence in the gap.55,56 In the present simulations, geometry optimizations and frequency calculations were performed with the periodic boundary condition by the Perdew−Burke− Ernzerhof57 functional as implemented in the Vienna ab initio simulation package.58 The semiempirical dispersion correction (D2)59 was considered because of the physisorption. Here the C6 and R0 for Au were set to 40.620 J nm6/mol and 1.772 Å, respectively.60 The parameters for other elements were taken from ref 59. The projector augmented-wave pseudopotentials,61,62 which took the scalar relativistic corrections into account, were used to represent the core electrons. Meanwhile, the wave functions were expanded by plane-wave basis sets with the energy cutoff of 400 eV. The Methfessel and Paxton method63 with a broadening factor of 0.1 eV was adopted for improving convergence of the electronic structure calculations, and the total energies were extrapolated to 0 K. The calculated lattice constant of Au at the Perdew−Burke− Ernzerhof functional level was 4.174 Å, which agrees well with the experimental value of 4.078 Å.64 With the calculated lattice constant, three layers of a (3√3 × 5) supercell of the Au(111) slab was used to simulate the Au surface. For the k-point sampling of the supercell, a 3 × 3 × 1 mesh of the Monkhorst− Pack grid65 was used, which generated five k-points in the irreducible Brillouin zone for the Brillouin zone integration. Enough vacuum layer (>35 Å) was added along the z axis in the supercell, and the dipole correction66,67 was used to avoid the artificial interactions. During the optimizations, the positions of the two bottom layers of the substrate were fixed to mimic the Au bulk, while all other degrees of freedom were allowed to relax. The force criterion was set to 0.02 eV/Å for optimizations. Furthermore, the central finite difference with displacement of 0.01 Å was used for frequency calculations with all atoms of the substrate fixed.

π 2cνs̃ 4Md 2|μ0 |2

r r r αpq = A pq + Bpq

(16)

where ΔErg is the corresponding vertical excitation energy. We should emphasize that the normal mode “k” is selected out because of the multimode nature of |vi⟩ and |vf⟩ and the integral ⟨vf |Qk|vi⟩. Besides, the sum-over-states of |vr⟩ is circumvented, which benefits practical calculations. Once μ0,k is obtained, the nonresonant Raman scattering intensity could be readily calculated by eq 12. It is noteworthy that the contributions of high orders with respect to Qk for both resonant and nonresonant processes can be obtained in a similar way by different truncations.24

(11)

2ϵ0

∂αeff 0 f Ei ⟨v |Q k|v i⟩ ∂Q k

with the matrix element of

If we consider the induced dipole as a classical oscillating dipole, the scattering intensity could be readily expressed as18,19,24 Is =

FP Mi

⟨v f |v r ⟩⟨v r |Q k|v i⟩

(15) 4988

DOI: 10.1021/acs.jctc.6b00592 J. Chem. Theory Comput. 2016, 12, 4986−4995

Article

Journal of Chemical Theory and Computation To calculate the Raman images, the coordinates of the adsorbate and all normal modes with frequency beyond 300 cm−1 in the periodic boundary condition calculations were extracted. All derivatives with respect to normal mode were calculated by the central finite difference68 with displacement of ΔQ = 1.0 × 10−3 amu1/2 bohr. Excited sates were calculated by the Gaussian suite of programs69 at the time-dependent density functional theory level with the B3LYP functional and Sadlej’s pVTZ basis set.70,71 Here effects of the weak interaction between adsorbate and substrate on excited states are neglected due to the limitations of the computational methods and, for nonresonant Raman processes, all singly excited states except core−hole excitations from oxygen and carbons were considered. Specifically, 2580 excited states were considered in our calculations for nonresonant Raman images. We extracted transition dipole moments and vertical excitation energies as well as the coefficients of excitation and deexcitation from calculations by the Gaussian program. The former two quantities are directly used for the calculations of α and αeff in eqs 8 and 17. For the electronic parts in eqs 8 and 17, the only unknown quantity is the value of ⟨ψg|q̂g(r)|ψr⟩, which contains the effects of SCP and thus is the key for the images. Using the coefficients of excitation (X) and de-excitation (Y), we have ⟨ψg |qĝ (r)|ψr ⟩ =

∑ (Xiar + Yiar)⟨ϕi|qĝ (r)|ϕa⟩

Ix = exp(αxtotGx 2 − αχK rχ , x 2 − ανL rν , x 2 − αx rD , x 2) ×

⟨ϕi|qĝ (r)|ϕa⟩ =

∑ cχicνa(∑ χν

KL

Gx =

∑ ∑ ∑ cαlmn,Dgαlmn,D

(24)

αχK rχ , x + ανL rν , x + αx rD , x αxtot

(25)

Here, all indices in the summation start out equal to zero and the prime (′) symbol in the summation represents that only the li + lj + lk = even terms were considered. Besides, the subscript of α highlights that, in eq 21, the exponent could have different values for different Cartesian components. By considering z as a Cartesian Gaussian function located at the original point with nz = 1 and α = 0, Iz becomes a four-center integral and can be calculated as Iz = exp(αztotGz 2 − αχK rχ , z 2 − ανL rν , z 2 − αz rD , z 2) ×

⎛ n K ⎞⎛ n L ⎞⎛ n ⎞⎛ 1 ⎞ L χ ⎟⎜ ν ⎟ n χK − ni (Gz − rν , z)nν − nj ⎟⎜ ⎟⎜ n ⎟⎜ ⎟(Gz − rχ , z) ⎝ ni ⎠⎝ nj ⎠⎝ k ⎠⎝ nl ⎠

∑ ′⎜⎜ ni nj

× (Gz − rD , z)n − nk (Gz)1 − nl

(ni + nj + nk + nl − 1)!! (2αztot)(ni + nj + nk + nl)/2

where αtot z and Gz could be calculated from eq 25 by replacing x with z. Alternatively, the more efficient recurrence formula would be adopted for the three-center and four-center integrals. Specifically, for the three-center integrals, the recurrence formula reads72

(19)

la ⟨la − 1|lD|lb⟩ 2αxtot l l + btot ⟨la|lD|lb − 1⟩ + Dtot ⟨la|lD − 1|lb⟩ 2αx 2αx

⟨la + 1|lD|lb⟩ = (Gx − rχ , x)⟨la|lD|lb⟩ +

(27)

where the integral over s functions is

lmn where gα,D is a Gaussian function localized at the center rD with exponent α, which can be written as

⟨0a , x|0D , x|0b , x⟩ = exp(αxtotGx 2 − αχK rχ , x 2 − ανL rν , x 2 − αx rD , x 2)

2

gαlmn = (x − xD)l (y − yD )m (z − zD)n e−α(r − rD) ,D

π αxtot

(28)

(21)

For the four-center integrals, we have the recurrence formula

and clmn α,D is the corresponding coefficient. We should emphasize that, here, rD is not necessarily equal to the tip position and all parameters in eq 20 could be obtained by fitting the realistic electric field distribution. Therefore, the integral

na ⟨na − 1|znD|nb⟩ 2αztot n n + btot ⟨na|znD|nb − 1⟩ + Dtot ⟨na|z(nD − 1)|nb⟩ 2αz 2αz 1 + ⟨na|nD|nb⟩ 2αztot (29)

⟨na + 1|znD|nb⟩ = (Gz − r χ , z)⟨na|znD|nb⟩ +

(22)

has to be estimated. In Raman images, we only considered the zz components in eqs 8 and 17. Thus, the integral eq 22 could be factored as I = IxIyIz

(26)

(20)

D l ,m,n α

I = ⟨g χK |qĝ αlmn |g L ⟩ ,D ν

π αztot

nk nl

where cχi and cνa are molecular orbital coefficients, gKχ and gLν are primitive Cartesian Gaussian functions, and DKχ and DLν are corresponding contraction coefficients. For the sake of computations, we expand arbitrary amplitude distribution by a set of Gaussian functions as g (r) =

(2αxtot)(li + l j+ lk)/2

αxtot = αχK + ανL + αx

where ϕi and ϕa represent occupied and unoccupied molecular orbitals, respectively. In eq 18, the integral could be calculated prior. In the representation of atomic orbital basis sets, it can be calculated as DχK DνL⟨g χK |qĝ (r)|gνL⟩)

(li + l j + lk − 1)!!

where rχ (rν) is the center of the atomic basis set χ (ν) and

(18)

ia

⎛l K ⎞⎛l L ⎞⎛ l ⎞ K L ⎟(Gx − rχ , x)l χ − li (Gx − rν , x)lν − l j (Gx − rD , x)l − lk l l ⎝ ⎠ k l ⎝ i ⎠⎝ j ⎠

∑ ′⎜⎜ χ ⎟⎟⎜⎜ ν ⎟⎟⎜ lil jlk

×

π αxtot

where the integral over s functions is ⟨0a , z|z 0D , z|0b , z⟩ = Gz⟨0a , z|0D , z|0b , z⟩

(30)

Again, the integral ⟨0a,z|0D,z|0b,z⟩ could be calculated from eq 28 by replacing x with z. In practical calculations, the horizontal recurrence relation73,74 was implemented to accelerate the computations. In addition, the transformation between

(23)

Here Ix and Iy are three-center overlap integrals that could be readily calculated as72 (hereafter taking Ix as an example) 4989

DOI: 10.1021/acs.jctc.6b00592 J. Chem. Theory Comput. 2016, 12, 4986−4995

Article

Journal of Chemical Theory and Computation Cartesian and pure spherical harmonic Gaussian functions75 was also performed if necessary. It should be noted that, for the same integrals, other recurrence relations based on the Hermite Gaussian76 or the hybrid scheme77 could be easily extended. For the vibrational parts in eq 8, the FC and HT integrals were evaluated with a widely used model for the excited state potential energy surface called linear coupling model (LCM).78 In this model the mode mixing effect (Duschinsky rotation effect) and frequency changes between the ground state and excited states were neglected which, thus, greatly simplified the quantum mechanical calculations since no geometrical optimization and vibrational analysis on the excited states were required. In LCM, the multimode FC and HT integrals can be computed by N independent single mode integrals as78 ⟨ν r|ν i⟩ =

is consistent with the distance between the tip and the substrates in experiments.15,81 The z position of the scanning plane was set to around 2 Å above the highest position of the adsorbate. To illustrate the resonant Raman processes, the light energy matched the first vertical excitation energy of 5.98 eV (207 nm) and the first 10 excited states were considered. Meanwhile, Γ was set to 100 cm−1 and only fundamental vibrations were considered. Our tested calculations show that 10 excited states are adequate for the resonant conditions. Besides, for nonresonant Raman spectra, the incident light was set to 632.8 nm. For images, we used a 81 × 81 grid with 0.2 Å step for the scanning plane, which generated 6561 individual points. For clarity, we used the relative intensity of different grid points to construct the Raman images. In other words, only the contrast was given for theoretical results; such treatment is the same as those reported in experiments.15 In addtion, the scale facors with respect to the maximum intensity of the CH stretching mode were also labeled in each image. All Raman images were calculated by the “first-principles approaches for surface- and tip-enhanced Raman scattering (FASTERS)” program.82 In addition, calculated spectra were convoluted by the Lorentzian function with fwhm of 10 cm−1. We should stress that, here, the resonant Raman images are only for conceptual demonstration, because it is difficult to generate plasmon in the deep-ultraviolet range. On the other hand, the nonresonant Raman images are more realistic.

∏ ⟨νkr|νki⟩ k

r

i

⟨ν |Q k|ν ⟩ = ⟨νkr|Q k|νki⟩ ∏ ⟨νlr|νli⟩ (31)

l≠k

νrk ≤ 79

In the case of be calculated as

i

νik,

the single mode FC integral

r

⟨νkr|νki⟩ = (− 1)νk − νk e−ωkΔQ k ⎞ ⎛ ωk ×⎜ ΔQ k ⎟ ⎠ ⎝ 2ℏ

2

⟨νrk|νik⟩

can

/4ℏ

νki − νkr

νr νkr! k νki! j = 0



⎛ ω ΔQ 2 ⎞ j νki! ⎜− k 2ℏ k ⎟ ⎝ ⎠ (νkr − j)! (νki − νkr + j)! j!

(32)

4. RESULTS AND DISCUSSION The optimized structure of AC adsorbed on the Au(111) surface was depicted in Figure 1a. We have noted that the

Here ωk and ΔQk are the frequency of mode k and displacement of the potential energy surfaces along normal mode Qk, respectively. ΔQk can be computed as ΔQ k = −

1 ∂ΔErg ωk 2 ∂Q k

(33)

where the derivation can be obtained from either the force differences between ground and excited states or the central finite difference method. When νrk > νik, ⟨νrk|νik⟩ can be calculated from the relationship79 r

i

⟨νkr|νki⟩ = ( −1)νk − νk ⟨νki|νkr ⟩

(34)

The single mode HT integral ⟨νrk|Qk|νik⟩ can be obtained by expressing Qk in terms of the lowering (a) and raising (a†) operators as78 Qk =

Then

ℏ (a + a†) 2ωk

⟨νrk|Qk|νik⟩

⟨νkr|Q k|νki⟩ =

(35)

can be rewritten as ℏ ( νki ⟨νkr|νki − 1⟩ + 2ωk

νki + 1 ⟨νkr|νki + 1⟩) (36)

which can then be obtained from the single mode FC integrals. The vibrational integrals were calculated with the DYNAVIB software.80 In practical simulations, we considered one s-type Gaussian function for g. It is noted that the theoretical limit of plasmonic confinement is determined by the Thomas−Fermi screening length, which can be down to 1 Å.17 Due to the small size of the adsorbate, unless explicitly specified, we set the fwhm for the x and y components of g to the theoretical limit of 1 Å.17,41 Meanwhile, the fwhm of the z component was set to 5 Å, which

Figure 1. (a) Sideview (left) and topview (right) for the optimized structure of an AC adsorbed on the Au(111) surface. Yellow, cyan, and white balls represent Au, C, and H atoms, respectively. PR, PC, and PL represent three tip positions. (b) Calculated resonant (top, red line) and nonresonant (bottom, black line) Raman spectra of AC with uniform EM field. (c) Representative vibrational modes for vCC, vCH, and vOH, respectively. The values are the corresponding frequencies in cm−1. 4990

DOI: 10.1021/acs.jctc.6b00592 J. Chem. Theory Comput. 2016, 12, 4986−4995

Article

Journal of Chemical Theory and Computation

Figure 2. (a) Calculated resonant Raman spectra of an AC adsorbed on the Au(111) surface when the center of SCP locates at PR, PC, and PL (from top to bottom), respectively. The values represent the scale factor for the corresponding spectrum. (b) Calculated resonant Raman images for vCC, vCH, and vOH (from left to right), respectively. The solid lines represent the skeleton of AC, and the values represent the scale factor for the corresponding image.

Figure 4. Calculated nonresonant linear (top) and nonlinear (bottom) Raman images of an AC adsorbed on the Au(111) surface for vCC, vCH, and vOH (from left to right), respectively. The x and y components of fwhm are set to (a) 5 and (b) 20 Å. The solid lines represent the skeleton of AC, and the values represent the scale factor for the corresponding image.

the substrate with a angle of around 37° with respect to the surface normal. We first report the resonant and nonresonant Raman spectra with the uniform EM (shown in Figure 1b), which would be good references for Raman images. One could immediately notice the differences between the resonant and nonresonant spectra. For instance, the bands higher than 2600 cm−1 are suppressed in the resonant spectrum. On the other hand, the CC stretching mode (vCC) located at 1566 cm−1 is significant in both resonant and nonresonant cases. Thus, we consider this mode for the simulation of Raman images. All the modes related to CH2 are located in the frequency region of 2600− 3200 cm−1. In the nonresonant spectrum, the four significant bands in this region correspond to the symmetrical C−H stretching, asymmetrical C−H stretching of the CH2 connected to the OH group, and symmetrical and asymmetrical C−H stretching of the CH2 connected to the vinyl group from low to high frequency. Here, we select the symmetrical stretching mode of the CH2 connected to the vinyl group (vCH) that locates at 2965 cm−1 to be the representative C−H stretching modes. It should be stressed that this mode has major component of C−H against the substrate as shown in Figure 1c. In the simulation of Raman images, the O−H stretching mode (vOH) located at 3548 cm−1 is also considered. We also noted that FC and HT terms both contribute to the resonant spectrum because of the small electronic transitions in the z direction. We then move on to the Raman spectra with SCP and the Raman images. For the spectra, we considered three representative positions for the center of the SCP, namely, PL, PC, and PR. These positions are above the left carbon atom of the vinyl group, the center carbon of the CH2 group, and the right oxygen atom of the OH group, respectively (see Figure 1a for details).

Figure 3. (a) Calculated nonresonant Raman spectra of an AC adsorbed on the Au(111) surface when the center of SCP locates at PR, PC, and PL (from top to bottom), respectively. The values represent the scale factor for corresponding spectrum. (b) Calculated nonresonant Raman images for vCC, vCH, and vOH (from left to right), respectively. The solid lines represent the skeleton of AC, and the values represent the scale factor for the corresponding image.

isolate AC has a weak intramolecular OH···π interaction.83 However, we considered the configuration that has stronger interaction with the substrate, which is in good agreement with previous studies for AC adsorbed on the same substrate.54 Specifically, the CC double bond adsorbed on the quasi-top site with the center around 2.6 Å above the substrate and a tilting angle of around 11°. The other two carbon atoms and the oxygen atom have almost the same z position (around 3.1 Å above the substrate). Meanwhile, the O−H bond points toward 4991

DOI: 10.1021/acs.jctc.6b00592 J. Chem. Theory Comput. 2016, 12, 4986−4995

Article

Journal of Chemical Theory and Computation

where Ms is the enhancement factor of the probe SCP, ωs is the resonant frequency, and the duration is 2 2 ln 2 τs . Here we assume Es has the same spatial amplitude distribution of E as shown in eq 1. As a result, the amplitude of induced resonant nonlinear dipole in the “two-state model” could be calculated as

The calculated resonant Raman spectra as well as the Raman images of the three representative vibrational modes (shown in Figure 1c) were depicted in Figure 2. It could be observed that although the spectra from PR and PL are both similar to the spectrum with uniform EM field, the resonant Raman spectra can be dependent on the positions of the SCP center. For instance, the high frequency modes are promoted in the spectrum from PC. This should be attributed to the fact that the HT term has significant contributions to the total intensity as shown in the conventional resonant spectrum. Surprisingly, the calculated images for different vibrational modes shown in Figure 2b all present a bright pattern located at the center CH2 group, although the modes of vCC and vOH are not directly associated with the center CH2. We also notice that the image of vCC has a bright pattern associated with the CC bond. We find that this result could be attributed to the fact that the first term in the HT term dominates the whole images, which is not frequency dependent as shown in eq 15. In this case, the resonant Raman images reflect the mixture of transition densities of the involving excited states. In order to obtain clear vibrational dependence in resonant Raman images, the second term of the HT term which is vibrationally dependent should be strong enough since both the FC term and the first term of the HT term are vibrationally independent. To find a system that has significant contribution from the second term of the HT term would be our future work. In contrast to the resonant conditions, with the inclusion of the SCP, the most important advantage of nonresonant Raman processes is vibrational selectivity. This feature can be immediately observed in the calculated nonresonant Raman spectra shown in Figure 3a. When the center of the SCP locates at PR, the vOH band is dominant in the whole spectra. On the other hand, when the center moves to PC, the vCH located at 2965 cm−1 is emerged with the suppression of the vOH band. Meanwhile, the other significant band located at 1376 cm−1 is the bending mode of ∠H−C−H in the center CH2 group. If the tip position locates at PL, the spectrum is dominated by the vCC band. These results manifest that, under nonresonant conditions, the individual vibrations can be selectively excited by changing the tip position, which has various potential applications in both fundamental and applied research. The selectivity is more clearly demonstrated in the calculated nonresonant Raman images as shown in Figure 3b. All images are located corresponding to their vibrational modes. It is nice to note that different shapes for individual images are presented in Figure 3b. Specifically, the image of vCH mode has quasiGaussian distribution, while the images of vCC and vOH both have stretching along the corresponding CC and O−H bonds. The nonspatial overlapping between different modes is consistent with the calculated spectra, which is also the prerequisite of the vibrational selection. In experiments, the generation of larger SCP is more feasible. However, for small adsorbates, the resolution would be decreased in this condition. One may expect the nonlinear effects could improve the resolution of the images. In fact, the nonlinear effects have significant contribution as reported in previous experimental Raman images.15 In the experiments, a broadband probe SCP is involved. The electric field of such broadband probe SCP could be expressed as

μ0NL = μ0,NL + μ0,NL A B where = μ0,NL A

Mi 2Ms|Ei0|3 2π τs ℏ3 ×

Es(r, t ) =

MsEi0 g(r) e−t /2τs

2

e

−ı ωst

⟨ψr|z|̂ ψg ⟩⟨ψr|zĝ (r)|ψg ⟩3

∑ ⟨v f |vr⟩⟨vr|vi⟩⟨vi|vr′⟩⟨vr′|v f ⟩λ(3) + NRT vr , vr′

(39)

= μ0,NL B

Mi 2 Ms 3

ℏ ×

|Ei0|3 2π τs[

∂⟨ψr|z|̂ ψg ⟩ ∂Q k

⟨ψr|zĝ (r)|ψg ⟩3

∑ ⟨v f |Q k|v r⟩⟨vr|vi⟩⟨vi|vr′⟩⟨vr′|v f ⟩λ(3)

vr , vr′

+ ⟨ψr|r |̂ ψg ⟩ ×

∂⟨ψr|zĝ (r)|ψg ⟩ ∂Q k

⟨ψr|zĝ (r)|ψg ⟩2

∑ (⟨v f |v r⟩⟨v r|Q k|vi⟩⟨vi|vr′⟩⟨vr′|v f ⟩

vr , vr′

+ ⟨v f |v r ⟩⟨v r|v i⟩⟨v i|Q k|v r ′⟩⟨v r ′|v f ⟩ + ⟨v f |v r ⟩⟨v r|v i⟩⟨v i|v r ′⟩⟨v r ′|Q k|v f ⟩)λ(3)] + NRT (40)

Here only the z component is considered and the definition

where ωo is the frequency of the detecting signal. Note that γ is the damping factor for different vibrational states in the same electronic state and three processes including stimulated Raman84,85 as well as two hot luminescence terms (hot luminescences I and II86) are considered in the definition of λ(3). We refer to ref 26 for the step-by-step derivation of eq 38. To estimate eq 38, we have to calculate two vibrational indices vr and vr′ in the summation, which is time-consuming. Thus, we do not calculate the nonlinear resonant Raman images here. On the other hand, because of the orthogonal of the vibrational integral, the amplitude of induced nonresonant nonlinear dipole could be calculated as26 μ0,NL =λ k

′ f ∂αeff ∂αeff ⟨v |Q k|v i⟩2 ∂Q k ∂Q k

(42)

where ′ pq = 2 ∑ αeff, r′

2

(38)

⟨ψg |pĝ (r)|ψr ′⟩⟨ψr ′|qĝ (r)|ψg ⟩ ΔEr ′ g

(43)

and λ is independent of SCP position

(37) 4992

DOI: 10.1021/acs.jctc.6b00592 J. Chem. Theory Comput. 2016, 12, 4986−4995

Article

Journal of Chemical Theory and Computation 2 2

M 2M |E 0|3 2π τs e−(ω − ωs − ωk −ı γ ) τs /2 λ= i s i ℏ ωo − ω + ωk + iγ

(5) Stöckle, R. M.; Suh, Y. D.; Deckert, V.; Zenobi, R. Nanoscale Chemical Analysis by Tip-Enhanced Raman Spectroscopy. Chem. Phys. Lett. 2000, 318, 131−136. (6) Anderson, M. S. Locally Enhanced Raman Spectroscopy with an Atomic Force Microscope. Appl. Phys. Lett. 2000, 76, 3130−3132. (7) Hayazawa, N.; Inouye, Y.; Sekkat, Z.; Kawata, S. Metallized Tip Amplification of Near-Field Raman Scattering. Opt. Commun. 2000, 183, 333−336. (8) Pettinger, B.; Pucardi, G.; Schuster, R.; Ertl, G. Surface Enhanced Raman Spectroscopy: Towards Single Molecule Spectroscopy. Electrochemistry (Tokyo, Jpn.) 2000, 68, 942−949. (9) Hartschuh, A.; Sánchez, E. J.; Xie, X. S.; Novotny, L. HighResolution Near-Field Raman Microscopy of Single-Walled Carbon Nanotubes. Phys. Rev. Lett. 2003, 90, 095503. (10) Anderson, N.; Hartschuh, A.; Cronin, S.; Novotny, L. Nanoscale Vibrational Analysis of Single-Walled Carbon Nanotubes. J. Am. Chem. Soc. 2005, 127, 2533−2537. (11) Steidtner, J.; Pettinger, B. Tip-Enhanced Raman Spectroscopy and Microscopy on Single Dye Molecules with 15 nm Resolution. Phys. Rev. Lett. 2008, 100, 236101. (12) Ichimura, T.; Fujii, S.; Verma, P.; Yano, T.; Inouye, Y.; Kawata, S. Subnanometric Near-Field Raman Investigation in the Vicinity of a Metallic Nanostructure. Phys. Rev. Lett. 2009, 102, 186101. (13) Yano, T.-a.; Verma, P.; Saito, Y.; Ichimura, T.; Kawata, S. Pressure-Assisted Tip-Enhanced Raman Imaging at a Resolution of a Few Nanometres. Nat. Photonics 2009, 3, 473−477. (14) Stadler, J.; Schmid, T.; Zenobi, R. Nanoscale Chemical Imaging Using Top-Illumination Tip-Enhanced Raman Spectroscopy. Nano Lett. 2010, 10, 4514−4520. (15) Zhang, R.; Zhang, Y.; Dong, Z. C.; Jiang, S.; Zhang, C.; Chen, L. G.; Zhang, L.; Liao, Y.; Aizpurua, J.; Luo, Y.; Yang, J. L.; Hou, J. G. Chemical Mapping of a Single Molecule by Plasmon-Enhanced Raman Scattering. Nature 2013, 498, 82−86. (16) Chen, C.; Hayazawa, N.; Kawata, S. A 1.7 nm Resolution Chemical Analysis of Carbon Nanotubes by Tip-Enhanced Raman Imaging in the Ambient. Nat. Commun. 2014, 5, 3312. (17) Atkin, J. M.; Raschke, M. B. Techniques: Optical Spectroscopy Goes Intramolecular. Nature 2013, 498, 44−45. (18) Albrecht, A. C. On the Theory of Raman Intensities. J. Chem. Phys. 1961, 34, 1476−1484. (19) Neugebauer, J.; Reiher, M.; Kind, C.; Hess, B. A. Quantum Chemical Calculation of Vibrational Spectra of Large MoleculesRaman and IR Spectra for Buckminsterfullerene. J. Comput. Chem. 2002, 23, 895−910. (20) Masiello, D. J.; Schatz, G. C. Many-Body Theory of SurfaceEnhanced Raman Scattering. Phys. Rev. A: At., Mol., Opt. Phys. 2008, 78, 042505. (21) Lee, S.-Y.; Zhang, D.; McCamant, D. W.; Kukura, P.; Mathies, R. A. Theory of Femtosecond Stimulated Raman Spectroscopy. J. Chem. Phys. 2004, 121, 3632−3642. (22) Sakurai, J. J. Advanced Quantum Mechanics; Addison-Wesley: Reading, MA, USA, 1967; pp 1−334. (23) Scully, M. O.; Zubairy, M. S. Quantum Optics, 1st ed.; Cambridge University Press: Cambridge, U.K., 1997; pp 1−630. (24) Long, D. A. The Raman Effect: A Unified Treatment of the Theory of Raman Scattering by Molecules; Wiley: Chichester, U.K., 2002; pp 1− 584, DOI: 10.1002/0470845767. (25) Le Ru, E.; Etchegoin, P. Principles of Surface-Enhanced Raman Spectroscopy: and Related Plasmonic Effects, 1st ed.; Elsevier: Amsterdam, 2009; pp 1−663. (26) Duan, S.; Tian, G.; Ji, Y.; Shao, J.; Dong, Z. C.; Luo, Y. Theoretical Modeling of Plasmon-Enhanced Raman Images of a Single Molecule with Subnanometer Resolution. J. Am. Chem. Soc. 2015, 137, 9515−9518. (27) Sun, M.; Zhang, Z.; Chen, L.; Sheng, S.; Xu, H. Plasmonic Gradient Effects on High Vacuum Tip-Enhanced Raman Spectroscopy. Adv. Opt. Mater. 2014, 2, 74−80.

(44)

Here we only consider the fundamental vibrations and the stimulated Raman process. Again, the zero-frequency limit is used. The calculated nonresonant linear and nonlinear Raman images with a larger SCP (5 and 20 Å) are depicted in Figure 4. As expected, the larger SCP gives more defocused images, although nonlinear effects indeed improve the resolution. Besides, we also find that the nonlinear images cannot give more molecular details than their linear counterparts. For instance, the feature of images for the CC and OH vibrations disappears. In other words, all images have quasi-Gaussian shape as shown in Figure 4.

5. CONCLUSION We have presented a general theory for modeling both resonant and nonresonant Raman images under the localized plasmonic field. The efficient recursive algorithm for required integrals is also derived. The resonant and nonresonant Raman images for the AC molecule adsorbed on the Au surface are comprehensively investigated with the newly developed method. The different contributions of FC and HT terms to the Raman images are clearly identified for both resonant and nonresonant cases. Our developed theoretical and computational methods can be easily interfaced with existing quantum chemical codes, which can lead to broad applications in simulations of Raman images of different systems.



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Author Contributions ∥

S.D. and G.T. contributed equally to this work.

Funding

This work was supported by the National Natural Science Foundation of China (Grant No. 21421063), the “Strategic Priority Research Program” of the Chinese Academy of Sciences (Grant No. XDB01020200), and the Swedish Research Council (VR). Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS The Swedish National Infrastructure for Computing (SNIC) is acknowledged for computer time.



REFERENCES

(1) Chalfie, M.; Tu, Y.; Euskirchen, G.; Ward, W. W.; Prasher, D. C. Green Fluorescent Protein as a Marker for Gene Expression. Science 1994, 263, 802−805. (2) Hell, S. W.; Wichmann, J. Breaking the Diffraction Resolution Limit by Stimulated Emission: Stimulated-Emission-Depletion Fluorescence Microscopy. Opt. Lett. 1994, 19, 780−782. (3) Betzig, E.; Patterson, G. H.; Sougrat, R.; Lindwasser, O. W.; Olenych, S.; Bonifacino, J. S.; Davidson, M. W.; Lippincott-Schwartz, J.; Hess, H. F. Imaging Intracellular Fluorescent Proteins at Nanometer Resolution. Science 2006, 313, 1642−1645. (4) Rust, M. J.; Bates, M.; Zhuang, X. Sub-Diffraction-Limit Imaging by Stochastic Optical Reconstruction Microscopy (STORM). Nat. Methods 2006, 3, 793−796. 4993

DOI: 10.1021/acs.jctc.6b00592 J. Chem. Theory Comput. 2016, 12, 4986−4995

Article

Journal of Chemical Theory and Computation (28) Meng, L.; Yang, Z.; Chen, J.; Sun, M. Effect of Electric Field Gradient on Sub-nanometer Spatial Resolution of Tip-Enhanced Raman Spectroscopy. Sci. Rep. 2015, 5, 9240. (29) Zhang, C.; Chen, B.-Q.; Li, Z.-Y. Optical Origin of Subnanometer Resolution in Tip-Enhanced Raman Mapping. J. Phys. Chem. C 2015, 119, 11858−11871. (30) Jiang, S.; Zhang, Y.; Zhang, R.; Hu, C.; Liao, M.; Luo, Y.; Yang, J.; Dong, Z.; Hou, J. G. Distinguishing Adjacent Molecules on a Surface Using Plasmon-Enhanced Raman Scattering. Nat. Nanotechnol. 2015, 10, 865−869. (31) Chen, X.; Moore, J. E.; Zekarias, M.; Jensen, L. Atomistic Electrodynamics Simulations of Bare and Ligand-Coated Nanoparticles in the Quantum Size Regime. Nat. Commun. 2015, 6, 8921. (32) Xu, H.; Wang, X.-H.; Persson, M. P.; Xu, H. Q.; Käll, M.; Johansson, P. Unified Treatment of Fluorescence and Raman Scattering Processes near Metal Surfaces. Phys. Rev. Lett. 2004, 93, 243002. (33) Lopata, K.; Neuhauser, D. Multiscale Maxwell-Schrödinger Modeling: A Split Field Finite-Difference Time-Domain Approach to Molecular Nanopolaritonics. J. Chem. Phys. 2009, 130, 104707. (34) Tian, G.; Liu, J.-C.; Luo, Y. Density-Matrix Approach for the Electroluminescence of Molecules in a Scanning Tunneling Microscope. Phys. Rev. Lett. 2011, 106, 177401. (35) Myers Kelley, A. Resonance Raman and Resonance HyperRaman Intensities: Structure and Dynamics of Molecular Excited States in Solution. J. Phys. Chem. A 2008, 112, 11975−11991. (36) Duan, S.; Tian, G.; Luo, Y. Visualization of Vibrational Modes in Real Space by Tip-Enhanced Non-Resonant Raman Spectroscopy. Angew. Chem., Int. Ed. 2016, 55, 1041−1045. (37) Ma, H.; Liu, J.; Liang, W. Time-Dependent Approach to Resonance Raman Spectra Including Duschinsky Rotation and Herzberg-Teller Effects: Formalism and Its Realistic Applications. J. Chem. Theory Comput. 2012, 8, 4474−4482. (38) Egidi, F.; Bloino, J.; Cappelli, C.; Barone, V. A Robust and Effective Time-Independent Route to the Calculation of Resonance Raman Spectra of Large Molecules in Condensed Phases with the Inclusion of Duschinsky, Herzberg-Teller, Anharmonic, and Environmental Effects. J. Chem. Theory Comput. 2014, 10, 346−363. (39) Jackson, J. D. Classical Electrodynamics, 3rd ed.; Wiley: New York, 1998; pp 1−808. (40) Zhu, W.; Esteban, R.; Borisov, A. G.; Baumberg, J. J.; Nordlander, P.; Lezec, H. J.; Aizpurua, J.; Crozier, K. B. Quantum Mechanical Effects in Plasmonic Structures with Subnanometre Gaps. Nat. Commun. 2016, 7, 11495. (41) Kong, J. A. Electromagnetic Wave Theory, 2nd ed.; Wiley: New York, 1990; pp 1−718. (42) Archambault, A.; Marquier, F.; Greffet, J.-J.; Arnold, C. Quantum Theory of Spontaneous and Stimulated Emission of Surface Plasmons. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 82, 035411. (43) Shankar, R. Principles of Quantum Mechanics, 2nd ed.; Plenum Press: New York, 1994; pp 1−676. (44) Zettili, N. Quantum Mechanics: Concepts and Applications, 1st ed.; Wiley: Chichester, U.K., 2001; pp 1−649. (45) Blume, M. Magnetic Scattering of X Rays (Invited). J. Appl. Phys. 1985, 57, 3615−3618. (46) Purcell, E. M. Spontaneous Emission Probabilities at Radio Frequencies. Phys. Rev. 1946, 69, 681. (47) Silverstein, D. W.; Jensen, L. Vibronic Coupling Simulations for Linear and Nonlinear Optical Processes: Theory. J. Chem. Phys. 2012, 136, 064111. (48) Baseggio, O.; Fronzoni, G.; Stener, M. A New Time Dependent Density Functional Algorithm for Large Systems and Plasmons in Metal Clusters. J. Chem. Phys. 2015, 143, 024106. (49) McHale, J. L. Molecular Spectroscopy, 1st ed.; Prentice Hall: Upper Saddle River, NJ, USA, 1998; pp 1−463. (50) Ringler, M.; Schwemer, A.; Wunderlich, M.; Nichtl, A.; Kürzinger, K.; Klar, T. A.; Feldmann, J. Shaping Emission Spectra of Fluorescent Molecules with Single Plasmonic Nanoresonators. Phys. Rev. Lett. 2008, 100, 203002.

(51) Johansson, P.; Xu, H.; Käll, M. Surface-Enhanced Raman Scattering and Fluorescence near Metal Nanoparticles. Phys. Rev. B: Condens. Matter Mater. Phys. 2005, 72, 035427. (52) Jensen, L.; Aikens, C. M.; Schatz, G. C. Electronic Structure Methods for Studying Surface-Enhanced Raman Scattering. Chem. Soc. Rev. 2008, 37, 1061−1073. (53) Gong, Z.-Y.; Tian, G.; Duan, S.; Luo, Y. Significant Contributions of the Albrecht’s A Term to Nonresonant Raman Scattering Processes. J. Chem. Theory Comput. 2015, 11, 5385−5390. (54) Mullen, G. M.; Zhang, L.; Evans, E. J.; Yan, T.; Henkelman, G.; Mullins, C. B. Oxygen and Hydroxyl Species Induce Multiple Reaction Pathways for the Partial Oxidation of Allyl Alcohol on Gold. J. Am. Chem. Soc. 2014, 136, 6489−6498. (55) Yang, Z.; Li, Q.; Ruan, F.; Li, Z.; Ren, B.; Xu, H.; Tian, Z. FDTD for Plasmonics: Applications in Enhanced Raman Spectroscopy. Chin. Sci. Bull. 2010, 55, 2635−2642. (56) Kazemi-Zanjani, N.; Vedraine, S.; Lagugné-Labarthet, F. Localized Enhancement of Electric Field in Tip-Enhanced Raman Spectroscopy Using Radially and Linearly Polarized Light. Opt. Express 2013, 21, 25271−25276. (57) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (58) Kresse, G.; Furthmüller, J. Efficient Iterative Schemes for Ab Initio Total-Energy Calculations Using a Plane-Wave Basis Set. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 11169−11186. (59) Grimme, S. Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27, 1787−1799. (60) Amft, M.; Lebègue, S.; Eriksson, O.; Skorodumova, N. V. Adsorption of Cu, Ag, and Au Atoms on Graphene Including van der Waals Interactions. J. Phys.: Condens. Matter 2011, 23, 395001. (61) Blö chl, P. E.; Jepsen, O.; Andersen, O. K. Improved Tetrahedron Method for Brillouin-Zone Integrations. Phys. Rev. B: Condens. Matter Mater. Phys. 1994, 49, 16223−16233. (62) Kresse, G.; Joubert, D. From Ultrasoft Pseudopotentials to the Projector Augmented-Wave Method. Phys. Rev. B: Condens. Matter Mater. Phys. 1999, 59, 1758−1775. (63) Methfessel, M.; Paxton, A. T. High-Precision Sampling for Brillouin-Zone Integration in Metals. Phys. Rev. B: Condens. Matter Mater. Phys. 1989, 40, 3616−3621. (64) Lide, D. R., Ed. CRC Handbook of Chemistry and Physics, 87th ed.; Taylor and Francis: Boca Raton, FL, USA, 2007; pp 1-1-B-5. (65) Monkhorst, H. J.; Pack, J. D. Special Points for Brillouin-zone integrations. Phys. Rev. B 1976, 13, 5188−5192. (66) Neugebauer, J.; Scheffler, M. Adsorbate-Substrate and Adsorbate-Adsorbate Interactions of Na and K Adlayers on Al(111). Phys. Rev. B: Condens. Matter Mater. Phys. 1992, 46, 16067−16080. (67) Makov, G.; Payne, M. C. Periodic Boundary Conditions in Ab Initio Calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 1995, 51, 4014−4022. (68) Fornberg, B. Generation of Finite Difference Formulas on Arbitrarily Spaced Grids. Math. Comp. 1988, 51, 699−760. (69) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö ; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, Revision A.02; Gaussian: Wallingford, CT, USA, 2009. 4994

DOI: 10.1021/acs.jctc.6b00592 J. Chem. Theory Comput. 2016, 12, 4986−4995

Article

Journal of Chemical Theory and Computation (70) Bauernschmitt, R.; Ahlrichs, R. Treatment of Electronic Excitations within the Adiabatic Approximation of Time Dependent Density Functional Theory. Chem. Phys. Lett. 1996, 256, 454−464. (71) Sadlej, A. J. Medium-Size Polarized Basis Sets for High-Level Correlated Calculations of Molecular Electric Properties. Collect. Czech. Chem. Commun. 1988, 53, 1995−2016. (72) Obara, S.; Saika, A. Efficient Recursive Computation of Molecular Integrals over Cartesian Gaussian Functions. J. Chem. Phys. 1986, 84, 3963−3974. (73) Head-Gordon, M.; Pople, J. A. A Method for Two-Electron Gaussian Integral and Integral Derivative Evaluation Using Recurrence Relations. J. Chem. Phys. 1988, 89, 5777−5786. (74) Hamilton, T. P.; Schaefer, H. F. New Variations in TwoElectron Integral Evaluation in the Context of Direct SCF Procedures. Chem. Phys. 1991, 150, 163−171. (75) Schlegel, H. B.; Frisch, M. J. Transformation Between Cartesian and Pure Spherical Harmonic Gaussians. Int. J. Quantum Chem. 1995, 54, 83−87. (76) McMurchie, L. E.; Davidson, E. R. One- and Two-Electron Integrals Over Cartesian Gaussian Functions. J. Comput. Phys. 1978, 26, 218−231. (77) Gill, P. M. W.; Head-Gordon, M.; Pople, J. A. An Efficient Algorithm for the Generation of Two-Electron Repulsion Integrals Over Gaussian Basis Functions. Int. J. Quantum Chem. 1989, 36, 269− 280. (78) Macak, P.; Luo, Y.; Ågren, H. Simulations of Vibronic Profiles in Two-Photon Absorption. Chem. Phys. Lett. 2000, 330, 447−456. (79) Tian, G.; Luo, Y. Isomer-Dependent Franck-Condon Blockade in Weakly Coupled Bipyridine Molecular Junctions. J. Phys. Chem. C 2014, 118, 14853−14859. (80) Tian, G.; Duan, S.; Hua, W.; Luo, Y. DynaVib, Version 1.0; Royal Institute of Technology: Stockholm, Sweden 2012. (81) Qiu, X. H.; Nazin, G. V.; Ho, W. Vibrationally Resolved Fluorescence Excited with Submolecular Precision. Science 2003, 299, 542−546. (82) Duan, S.; Tian, G.; Luo, Y. First-Principles Approaches for Surface and Tip Enhanced Raman Scattering (FASTERS), Version 1.0; 2016; Royal Institute of Technology: Stockholm, Sweden. (83) Mackeprang, K.; Schrøder, S. D.; Kjaergaard, H. G. Weak Intramolecular OH···π Hydrogen Bonding in Methallyl- and Allylcarbinol. Chem. Phys. Lett. 2013, 582, 31−37. (84) Eckhardt, G.; Hellwarth, R. W.; McClung, F. J.; Schwarz, S. E.; Weiner, D.; Woodbury, E. J. Stimulated Raman Scattering From Organic Liquids. Phys. Rev. Lett. 1962, 9, 455−457. (85) Woodbury, E. J.; Ng, W. K. Ruby Laser Operation in the Near IR. Proc. IRE 1962, 50, 2347−2348. (86) Shen, Y. R. Distinction Between Resonance Raman Scattering and Hot Luminescence. Phys. Rev. B 1974, 9, 622−626.

4995

DOI: 10.1021/acs.jctc.6b00592 J. Chem. Theory Comput. 2016, 12, 4986−4995