Mixed atomistic–implicit quantum/classical approach to molecular nanoplasmonics Open Access

When metal nanoparticles (NPs) are irradiated with external radiation, coherent oscillations of conduction electrons, also named localized surface plasmons (LSPs), can be excited.^1–3 Most NP optical properties are related to LSP peculiar properties, which can be tuned by varying the NP’s size, shape, and chemical composition.⁴ Plasmonic nanostructures can enhance, control, or suppress properties of molecules interacting with light: these features are exploited in molecular nanoplasmonics.^5–10 A deep understanding of the phenomena that occur at the molecular and nanoscale in the presence of light can be achieved by exploiting multiscale hybrid techniques, which use different levels of description for molecules and plasmonic nanosystems.^5,11 With these methods, a reliable representation of both atomistic details and collective features, such as plasmons, in these complex systems can be achieved.

In principle, the ideal theoretical approach to molecular nanoplasmonics should rely on a quantum mechanical (QM) description of the whole system (molecule + nanostructured substrate). However, due to their unfavorable computational scaling, full ab initio methods are currently limited to small model systems (generally, $<$100 atoms), in contrast to the large size of plasmonic substrates, which typically comprise thousands or millions of atoms.^12–16 To address these limitations, hybrid multiscale QM/classical approaches have been developed,^11,17,18 where the molecular adsorbate is treated at the QM level, while the substrate’s plasmonic response is calculated by using classical electrodynamical methods, substantially reducing the computational cost. The most crude classical description of the plasmonic response of a metal NP is given by continuum implicit approaches, such as the Mie theory,¹⁹ the Finite Difference Time Domain (FDTD),²⁰ or the Boundary Element Method (BEM),^21–24 which completely disregard the NP atomistic nature. To overcome this limitation, fully atomistic approaches have been developed, such as Discrete Interaction Models (DIMs)^{17,18,25–30} and the Frequency-Dependent Fluctuating Charges and Fluctuating Dipoles (ωFQFμ) approach.^{11,24,31–39} ωFQFμ can correctly reproduce QM results for metal NPs even below the quantum limit (⁠$<$5 nm),³⁵ capture the plasmonic properties of systems featuring subnanometer junctions,³¹ defects,^34,36 treat bimetallic particles,³⁷ and colloidal nanostructures.³⁹ In addition, it describes noble metal NPs³⁵ and graphene-based structures under the same theoretical framework.³² Both families (implicit and atomistic) have pros and cons. Continuum models feature a favorable computational cost, which scales with the size of the NP surface, facilitating the description of large systems.^23,40–46 However, their implicit, non-atomistic nature fails at capturing the NP plasmonic response in specific configurations characterized by atomistic defects, sub-nanometer junctions, and sharp interfaces,^45,47 which are associated with huge enhancements of the electric field (the so-called hot-spots). Conversely, while atomistic approaches appropriately capture these features, they become computationally less efficient as the size of the system increases, because their computational cost scales with the number of atoms, i.e., the NP volume.³³ 

This work proposes a novel multiscale approach specifically designed to overcome the limitations associated with atomistic and continuum approaches. The method, which is here specified for noble metal NPs, describes the core with an implicit approach, by using the BEM method, while the NP surface is treated at the fully atomistic approach, by means of ωFQFμ. The resulting ωFQFμ-BEM approach constitutes, to the best of our knowledge, the first hybrid atomistic–continuum methodology to evaluate the optical response of plasmonic substrates within classical electrodynamics. Furthermore, ωFQFμ-BEM is coupled to a QM description of a molecular adsorbate described at the Time-Dependent Density Functional Theory (TDDFT) level, allowing for the calculation of molecular properties and signals in the vicinity of plasmonic nanostructures.^11,42,48

One of the most interesting aspects of molecular nanoplasmonics is the huge enhancement of the induced electric field in the NP surface proximity, which can drastically affect the electronic properties of molecular adsorbates. As a result, molecular spectral signals can be significantly enhanced, providing an invaluable platform for molecular sensing. The most diffuse technique that exploits this effect is Surface Enhanced Raman Scattering (SERS),^{2,6,10,14,49–57} where the molecular Raman signals are enhanced by several orders of magnitude, allowing single molecule detection. SERS is at present used in various applicative fields, such as catalysis, chemical biology, biophysics, and biomedicine.^58–65 For this reason, in this paper, QM/ωFQFμ-BEM is extended to compute SERS signals.

This paper is organized as follows: In Sec. II, ωFQFμ and BEM methods are briefly recalled, and the novel ωFQFμ-BEM and QM/ωFQFμ-BEM models are presented. After a brief section reporting on the computational details, the methods are validated by computing NPs optical properties and SERS spectra of pyridine with the novel approaches or employing reference fully atomistic methods (ωFQFμ and QM/ωFQFμ). Conclusions and future perspectives end the paper.

This section gives an overview of the theoretical background leading to the formulation of the QM/ωFQFμ-BEM approach. First, the integration of the ωFQFμ and BEM models is discussed, and the ωFQFμ-BEM equations are presented to replicate the plasmonic response of fully atomistic NPs as built by a continuum BEM core within an atomistic ωFQFμ shell (see Fig. 1). Then, the QM/ωFQFμ-BEM coupling is developed and specified to describe the SERS of a molecule, studied at the TDDFT level, placed in the vicinity of the plasmonic NP.

The Boundary Element Method (BEM) constitutes a reliable approach for solving the Polarizable Continuum Model equations applied to the study of plasmonic nanoparticles (PCM-NPs).⁴¹ There, the NP is treated as a homogeneous continuum dielectric by exploiting a classical electrodynamics formalism. Plasmonics is then governed by the frequency-dependent permittivity function chosen to describe the material and that of the media in which it is embedded.^{23,24,45,47,71} BEM involves the calculation of surface charges, which yield a global surface charge density that encapsulates the optical response of the system. BEM charges are calculated upon applying the model on a triangular-discretized mesh delimiting the NP surface. At each triangular centroid, a complex point charge that describes the optical response of the material is calculated. This enables the exploration of the plasmonic properties of the analyzed structures, which is influenced by the geometric features of the discretized mesh and by the material’s permittivity function.

In Sec. S1 of the supplementary material, further details on the F and S matrices are given, along with the description of each term in Eq. (10).

The idea of integrating ωFQFμ and BEM, giving rise to a model that is named ωFQFμ-BEM, arises from their close similarities. In particular, both formalisms are grounded in classical electrodynamics and rely on complex frequency-dependent quantities (charges/dipoles) to describe the optical response of metal NPs. However, these models present notable differences. First, the BEM response is rooted in the permittivity chosen to describe the material and lacks any atomistic description. Conversely, ωFQFμ is based on textbook concepts, such as the Drude model of conduction for describing intraband transitions in terms of fluctuating charges. Interband transitions are recovered through a set of fluctuating dipoles that are defined in terms of the experimental interband polarizability. By this, ωFQFμ can effectively decouple the physical processes influencing the optical response of the plasmonic substrate while also retaining its atomistic nature. However, ωFQFμ encounters scalability challenges, eventually leading to the exploration of alternative strategies to solve Eq. (8).^33,74 Remarkably, ωFQFμ-BEM can alleviate these limitations by representing NPs as built of an internal continuum BEM core and an external atomistic ωFQFμ shell, with substantial computational savings.

The terms in Eq. (14) are detailed in Sec. S1 of the supplementary material.

Once ωFQFμ-BEM has been formulated, it can be coupled to a QM description of a molecular system in a QM/Classical fashion,^{5,17,23,76–80} giving rise to the QM/ωFQFμ-BEM model. ωFQFμ-BEM describes the system’s response to external oscillating electric fields; therefore, it can be naturally translated into a linear response formalism. The coupling is done by following the same strategy that was exploited to formulate the two-layer QM/ωFQFμ approach.¹¹ 

In Eq. (16), X and Y are the excitation and de-excitation transition densities, respectively, while Γ is a phenomenological damping factor. ɛ is the molecular orbital energy, (ai|bj) denotes two-electron integrals, and the c_x, c_l parameters vary depending on the DFT functional employed. In addition, C^{QM/ωFQFμ-BEM} describes the polarization induced by the ωFQFμ-BEM charges and dipoles, which respond to the perturbed TDKS density.^11,17,84 On the right-hand-side, Q_ia = ⟨ϕ_i|V^α,pert(r, ω)|ϕ_a⟩ represents the expectation value of the perturbation operator polarized along α.

Ag structures are created by using the Atomic Simulation Environment (ASE) Python module v. 3.17.⁹⁰ In particular, Ag atoms are disposed in a Face-Centered Cubic (FCC) arrangement defined by a lattice parameter of 4.08 Å. Two morphologies are studied: spherical and icosahedral (Ih) NPs. For ωFQFμ calculations, the parameters defined in Ref. 35 are exploited. Detailed information regarding the geometrical characteristics of all atomistic structures is reported in Tables S2–S5 of the supplementary material.

The continuum BEM meshes representing the surfaces of both spherical and Ih NPs are constructed by using the GMSH software.⁹¹ The frequency-dependent permittivities proposed by Palik,⁹² Brendel-Bormann,⁹³ and Johnson and Christy⁹⁴ are exploited in BEM calculations.

ωFQFμ-BEM spherical and Ih structures are built as a core, described as an implicit spherical NP and treated at the BEM level, surrounded by an atomistic shell, treated at the ωFQFμ level. To validate the novel approach, we investigate the variation of the optical response of a spherical Ag nanostructure (radius 4 nm) as a function of the method parameters, namely, the core-to-shell distance, the thickness of the atomistic shell, the BEM tessellation, and the dielectric function to describe BEM response.

In particular, we study the optical response as a function of the minimum core-to-shell distance, which varies from 2.88 to 5.76 Å, corresponding to integer multiples of the nearest neighbor distance (2.88 Å) in the studied FCC lattices. For Ag Ih structures, the atomistic shell is defined, ensuring a minimum thickness of three atomic layers in the thinner region.

After validating ωFQFμ-BEM, we exploit the multiscale QM/ωFQFμ-BEM to calculate SERS signals of pyridine adsorbed on Ag Ih NPs (with radius varying from 1.9 to 3.8 nm). SERS spectra are computed by representing the Ag Ih NP at the fully atomistic ωFQFμ, implicit BEM, and hybrid ωFQFμ-BEM levels, thus allowing for validating the novel method by a robust comparison with state-of-the-art methods. In all cases, pyridine is placed longitudinally at a distance of 3 Å from the NP surface. Such distance is calculated from the nitrogen atom to the nearest NP atom/tessera centroid defining the tip of the atomistic/continuum surface of Ag Ih structures.

In all calculations, pyridine is described at the QM level utilizing the BP86 functional and a double-ζ-polarized DZP basis set, in agreement with previous studies.^11,86 TDDFT equations are solved by imposing a damping factor Γ = 0.01 eV.^11,86,95,96 $αk′$ and $γk′$ [see Eq. (19)] are calculated by a numerical differentiation scheme, with a constant step size of 0.001 Å.^11,97–99 Normal modes’ displacements of pyridine are evaluated in vacuo, without explicitly considering NP effects, which are expected to be small for the considered system.^11,86 The influence of the plasmon resonance on SERS signals is assessed by matching the incident light frequency for Raman with the Plasmon Resonance Frequency (PRF) of the NP, whose values are collected in Tables S5 and S6 in the supplementary material. Finally, SERS spectra are plotted using Lorentzian band-shapes characterized by a full width at half maximum (FWHM) of 4 cm⁻¹. QM/ωFQFμ-BEM, QM/ωFQFμ, and QM/BEM calculations are performed by using a locally modified version of the Amsterdam Modeling Suite (AMS).⁸¹ 

A summary of the parameters involved in ωFQFμ, BEM, ωFQFμ-BEM, and SERS calculations is given in Sec. S2.1 of the supplementary material.

In this section, the ωFQFμ-BEM and QM/ωFQFμ-BEM methods are validated by comparison with the reference ωFQFμ and QM/ωFQFμ approaches. In particular, ωFQFμ-BEM is first challenged to reproduce the plasmonic features of spherical NPs by analyzing how a variation of the model parameters affects the optical response. Then, QM/ωFQFμ-BEM is applied to compute SERS of pyridine adsorbed on complex-shaped Ih nanostructures.

In Fig. 2(a), we show a graphical representation of a spherical NP described at the ωFQFμ-BEM level, constructed by integrating a continuum BEM spherical core within an atomistic ωFQFμ spherical shell. ωFQFμ-BEM spherical NP geometries are defined by four parameters: the radius of the inner BEM core (r_BEM), the radii of the outer ωFQFμ shell (R), the BEM core-ωFQFμ shell distance (d), and the difference between the radii of the inner and outer shells (Δr). The computed optical response depends on the proper definition of such parameters, the dielectric function exploited to model the BEM portion, and the number of tesserae exploited to mesh the BEM region. It is worth noting that our approach is general, and the full BEM (r_BEM = R) or the full ωFQFμ (Δr = R) description can be easily recovered.

Let us first focus on the computed response for a specific geometry, i.e., a spherical NP with R = 40.00 Å, r_BEM = 20.26 Å, Δr = 13.98 Å, and d = 5.76 Å. 11 470 atoms constitute the atomistic shell, while the BEM core is defined by 2972 tesserae. The Brendel–Bormann frequency-dependent permittivity⁹³ is exploited to describe the BEM part.

Figure 2(b) compares normalized absorption spectra calculated for ωFQFμ-BEM and the reference, fully atomistic ωFQFμ (hereafter referred to as full-ωFQFμ). The normalization is performed with respect to the full-ωFQFμ PRF absorption. The full-ωFQFμ approach predicts a single absorption peak at 3.47 eV, whereas the ωFQFμ-BEM spectrum features two peaks centered at 2.50 and 3.45 eV. To analyze the physical origin of such a discrepancy, in Fig. 2(c), normalized absorption spectra of the selected spherical NP calculated at the ωFQFμ-BEM level are reported, together with the absorption spectrum of the BEM core and ωFQFμ shell. Spectra are normalized with respect to the maximum absorption of the ωFQFμ shell. The spherical shell spectrum shows a high-intensity peak centered at 3.24 eV and a low-intensity peak at 3.68 eV. These two peaks can be assigned to the bonding (|ω_s+⟩) and anti-bonding (|ω_s−⟩) plasmonic modes, respectively [see Fig. 2(d)], which are typical of plasmonic nanoshells.^24,100–102 A single plasmonic peak at about 3.43 eV (|ω_c⟩) is instead present in the spectrum of the BEM core. Such a band corresponds to the dipolar plasmonic excitation [see Fig. 2(d)].

The two peaks reported in the ωFQFμ-BEM spectrum of the spherical NP are thus due to the hybridization of isolated plasmonic modes from the core and shell, similar to what has been reported for other nanostructures.^102–106 More specifically, the ωFQFμ-BEM peak at 2.50 eV is associated with the hybrid |ω_s−⟩ − |ω_c⟩ mode, while the band centered at 3.45 eV corresponds to the |ω_s+⟩ + |ω_c⟩ mode combination [see Fig. 2(d)]. The presence of the low-energy peak is not unexpected because the approach discards any charge transfer between the two classical regions. This peak thus arises as an artifact of our model. However, it can be easily identified because it will always be the lowest energy peak computed in the spectrum (without the need for any benchmarking with a reference calculation), as predicted by hybridization theory.^102,105 In the following, the focus will thus shift to the |ω_s+⟩ + |ω_c⟩ band, whose PRF remarkably aligns with that predicted at the full-ωFQFμ level. For clarity, Fig. S1 of the supplementary material presents a comparison between the densities computed from the charges + dipoles distribution of the ωFQFμ-BEM |ω_s+⟩ + |ω_c⟩ band and the full-ωFQFμ PRF, illustrating a clear one-to-one correspondence.

As a final remark, a third |ω_s+⟩ − |ω_c⟩ mode arises from the hybridization of plasmonic modes. However, the small dipole moment associated with this mode makes it undetectable in the spectra.^102,105 Moreover, the relative intensities of the |ω_s−⟩ − |ω_c⟩ and |ω_s+⟩ + |ω_c⟩ modes are consistent with the expected hybridization behavior. In particular, the stronger intensity of the |ω_s+⟩ + |ω_c⟩ mode is expected due to the in-phase oscillation of the core and shell dipole moments.^102,105

This section first examines how a variation in the shell thickness and core–shell distance d influences the optical response of the metal NPs. In particular, the shell thickness Δr varies from 13.98 to 7.93 Å, with a constant step of ∼2.0 Å. Two core-to-shell distances are considered: d = 2.88 Å and d = 5.76 Å, corresponding to integer multiples of the nearest neighbor distance (2.88 Å) in the studied FCC lattice. In all cases, the BEM matrix is described by the Brendel–Bormann dielectric function.⁹³ 

Figures 3(a)–3(c) reports absorption cross sections of the spherical NP as a function of shell thickness Δr for d = 2.88 and d = 5.76 Å, respectively. Data are normalized with respect to the full-ωFQFμ PRF absorption. As a reference, the normalized absorption spectrum computed at the full-ωFQFμ level is reported in all panels. All ωFQFμ-BEM spectra feature an intense peak in the region 3.0–4.5 eV. Such a peak is associated with a dipolar plasmon, as commented in Sec. IV A 1. ωFQFμ-BEM computed PRFs remain almost constant as a function of the shell thickness and only present slight shifts of 0.02–0.04 eV with respect to ωFQFμ PRF (3.47 eV), regardless of the value of d. Differently, reducing the shell thickness systematically increases the absorption intensity mismatch at the PRF between ωFQFμ-BEM and the reference ωFQFμ data. Notably, the mismatch becomes more pronounced for structures with d = 2.88 Å. This behavior is probably due to numerical instabilities occurring when solving the ωFQFμ-BEM linear equation [see Eq. (14)] for the smaller d values, associated with the so-called polarization catastrophe raising for the increase of electrostatic/polarization interactions.¹⁰⁷ 

To further analyze the plasmonic response as a function of d, Figs. 3(b) and 3(d) graphically illustrate the differences in the computed charges + dipoles density at the ωFQFμ-BEM and full-ωFQFμ levels. In all cases, densities are computed at the PRF for d = 2.88 Å and d = 5.76 Å, respectively. This qualitative analysis permits us to evaluate the charges + dipoles density deviation near the NP surface of ωFQFμ-BEM with respect to the full-ωFQFμ atomistic structure, whose surface is graphically represented by green dots. The most remarkable differences are observed at the core region, and in particular at the BEM-ωFQFμ interface. This is expected and is due to the artificial boundary that is introduced in our multiscale method. However, since we are dealing with localized surface plasmons, we are particularly interested in a proper description of surface properties. Note that the differences at the NP surface are negligible for all methods, thus validating the novel methodology. The agreement between the two methods worsens by reducing the atomistic shell thickness, especially for d = 2.88 Å. This is again related to numerical instabilities in solving the ωFQFμ-BEM linear equation.

To conclude this section, we investigate the computational savings associated with the multiscale ωFQFμ-BEM as compared to a full ωFQFμ description of the spherical NP. To this end, the relative computational speed-up (in percentage) for single-frequency calculations (at the PRF of the systems) is reported in Table I as a function of d and Δr. We remark that all calculations exploit a constant number of tesserae to mesh the BEM core (∼2980 tesserae). The data reported in Table I clearly show a substantial computational saving, ranging from about 43% to about 77%, independent of d for a given value or Δr. This is expected because the number of atoms in the ωFQFμ shell decreases while the number of tesserae representing the BEM core remains constant. It is worth remarking that for the most accurate ωFQFμ-BEM partitioning (Δr = 13.98 Å, d = 5.76 Å), the speed-up of ωFQFμ-BEM is considerable, without losing accuracy as compared to a full ωFQFμ description. Therefore, the above model parameters represent the best compromise between accuracy and computational cost for ωFQFμ-BEM applications.

We now move to analyze the model parameters affecting the BEM core response, namely, the dielectric function and the number of tesserae. To this end, we consider three dielectric functions that are commonly exploited to describe Ag response, recovered from Brendel–Bormann (BB),⁹³ Palik,⁹² and Johnson and Christy (J&C).⁹⁴ Computed absorption cross sections for the selected spherical NP exploiting the three permittivity functions are normalized with respect to the full-ωFQFμ reference and graphically depicted in Figs. 4(a)–4(c), for d = 5.76 Å (see Fig. S3 of the supplementary material for additional values obtained with d = 2.88 Å). Absorption properties are reported as a function of the shell thickness Δr. The normalized full atomistic ωFQFμ absorption spectrum is also shown as a reference.

All computed spectra are characterized by a main plasmonic peak, which is associated with a dipolar localized surface plasmon. ωFQFμ-BEM computed PRFs (see also Table II) remain constant as a function of the shell thickness by exploiting Palik and J&C permittivity functions (3.48 eV), and remarkably the computed PRF is shifted by only 0.01 eV with respect to the reference ωFQFμ PRF (3.47 eV). As commented in Sec. IV A 2, by using the BB ɛ(ω), the PRF slightly shifts by decreasing the shell thickness (from 3.45 eV–Δr = 13.98 Å to 3.43 eV–Δr = 7.93 Å). However, also in this case, all data agree with the reference ωFQFμ PRF within the chemical accuracy (0.04 eV to 0.9 kcal/mol).

The absorption cross section at the PRF varies significantly depending on the chosen dielectric function. Notably, absorption intensities in very good agreement with the reference ωFQFμ model are obtained by employing BB and Palik dielectric functions. The ωFQFμ-BEM spectrum computed by using the Palik dielectric function presents a small shoulder at about 3.3 eV, which is also observed at the full BEM level (see Fig. S4 in the supplementary material) and is thus related to the use of this specific permittivity function. In contrast, a pronounced intensity mismatch as compared to the full-ωFQFμ absorption curve is observed if the J&C dielectric function is employed. Remarkably, such a discrepancy increases as the shell thickness decreases.

To quantitatively analyze surface-related near-field properties, ρ^error values as a function of the permittivity function and Δr are reported in Table II. Computed relative errors using BB and J&C permittivities consistently increase by reducing the atomistic shell thickness, unlike those obtained by exploiting the Palik dielectric function, which remains almost constant (about 9%). Notably, the best results (i.e., lowest errors) are obtained by employing BB, with the exception of the case Δr = 7.93 Å for which Palik gives the best results. J&C data substantially deviate from reference ωFQFμ values. On the other hand, BB emerges as the most robust permittivity function model to exploit in ωFQFμ-BEM applications.

As a final validation, we investigate the dependence of absorption spectra on the number of tesserae used to mesh the BEM core. In particular, we consider four different tessellations of the spherical BEM core (2990, 1948, 1018, and 390 tesserae), while keeping the other model parameters fixed to the best set resulting from the previous analysis (d = 5.76 Å, Δr = 13.98 Å, BB permittivity function). Numerical results are reported in Fig. S5 and Table S3 of the supplementary material. Interestingly, the discretization of the core only minimally affects the plasmonic response (spectra, density-differences plots, and ρ^error), while yielding substantial computational % speed-ups, up to 56.68% for 390 tesserae. To conclude, the best parameter set, which guarantees the best compromise between computational cost and accuracy, is defined by d = 5.76 Å, Δr = 13.98 Å, the BB permittivity function, and 390 tesserae.

As mentioned in Sec. I and specified in Sec. II B, ωFQFμ-BEM can be coupled to a QM Hamiltonian in a QM/Classical fashion and extended to compute SERS spectra. This section showcases the potentialities of the method to compute SERS spectra of pyridine near the tip of complex-shaped NPs. To this end, Ag icosahedral (Ih) NPs of increasing size (from 1.9 to 3.8 nm) are selected, and additional calculations modeling the molecule–NP system at the fully atomistic (QM/ωFQFμ) and fully implicit (QM/BEM) levels are discussed. A graphical representation highlighting the structural differences between the various methods is given in Fig. 5, panel (a), whereas geometrical parameters are shown in Table S4 in the supplementary material. To fully characterize ωFQFμ-BEM Ih systems and, similarly, spherical NPs, four parameters are exploited (the BEM radius r_BEM, the BEM-ωFQFμ distance d, the average thickness of the ωFQFμ shell Δr, and the NP radius R), while the equivalent full-ωFQFμ and full-BEM structures are described by a single radius parameter (R). Detailed data on the geometrical features and PRFs of the systems are gathered in Tables S4–S6 of the supplementary material.

In Fig. 5, the SERS signal of pyridine on Ag Ih NPs is studied as a function of the NP radius. Pyridine is adsorbed perpendicularly to the NP surface at a distance of 3 Å, with the N atom lying closest to the NP. In agreement with the preliminary analysis in Sec. IV A, the BB dielectric function is employed, and d is set to 5.76 Å. Figure 5(b) illustrates SERS intensities as a function of the NP radius, normalized with respect to the largest signal exhibited by each model. In all cases, this corresponds to the SERS signal of pyridine interacting with the largest Ih NP. Notably, for all models, only the signals associated with specific normal modes are enhanced (see Fig. S7 in the supplementary material for their graphical depiction). This is due to the fact that such vibrations are related to an orthogonal movement of pyridine atoms with respect to the NP surface. As a consequence, they experience the largest field gradient variation and are, therefore, preferentially enhanced, as has recently been shown by some of us.¹¹ The analysis of the results reveals significant similarities between the QM/ωFQFμ-BEM and QM/ωFQFμ approaches while also highlighting substantial discrepancies by using QM/BEM. In particular, the observed Raman peaks are consistent across all models when pyridine is adsorbed on the largest NPs (radius $>$2.7 nm). On the contrary, the smallest full-BEM structures (radius $<$2.5 nm) substantially deviate from the full atomistic picture [see Fig. 5(b)]. To better compare the three approaches, Fig. 5(c) presents normalized Raman spectra for the largest NP studied using each theoretical framework. The spectra depicted in Fig. 5(c) show an almost perfect alignment of the relative Raman intensities between ωFQFμ-BEM and full-ωFQFμ NPs. In contrast, substantial discrepancies are evident for full-BEM structures, particularly for the Raman peaks centered at 1024, 1207, 1460, and 1568 cm⁻¹.

To comprehensively analyze the performance of the various models, Fig. 5(d) collects AEF and MEF descriptors calculated as a function of the NP radius. For all methods, such indices are monotonically increasing. The reference QM/ωFQFμ approach predicts a maximum AEF and MEF of about 5000 and 12 000, respectively, thus displaying an overall enhancement of about 10³. QM/BEM shows a substantial quantitative mismatch, providing AEF and MEF that are completely overestimated by several orders of magnitude (about 10⁸). Conversely, QM/ωFQFμ-BEM presents a clear trend of quantitative similarity with the reference fully atomistic results, predicts AEF and MEF of the same order of magnitude, and only deviates from the reference by a factor of about 1.35 for the largest structure. In addition, in Fig. S8 of the supplementary material, normalized gas-phase Raman and SERS spectra computed at PRF are reported for the largest nanoparticle studied using each formalism.

It is worth remarking that the QM/BEM trends arising from Fig. 5 can be substantially affected by varying the dielectric functions (J&C and Palik). These results are reported in Sec. S2.3 of the supplementary material and show that J&C permittivity yields a complete modification of the spectral response for the largest NP (see Figs. S9 and S12 of the supplementary material), also affecting computed AEF and MEF values, further magnifying the quantitative discrepancies compared to the fully atomistic ωFQFμ description. Conversely, the QM/ωFQFμ-BEM spectra remain consistent regardless of the chosen parameters but are crucially affected by a variation in the BEM-ωFQFμ distance d. In fact, by setting d = 2.88 Å, inconsistent trends in both AEF and MEF are reported as the system size increases (see Figs. S11–S13 of the supplementary material). This behavior completely aligns with our findings from the study on ρ^error of spherical ωFQFμ-BEM NPs (see Sec. IV A).

To remove any arbitrariness in selecting the minimum NP-pyridine distance for the full-BEM systems, we have performed distance-dependent studies of AEF and MEF for pyridine adsorbed on the tip of icosahedral nanoparticles with a 3 nm radius. These calculations were performed at the full-BEM, full-ωFQFμ, and ωFQFμ-BEM (characterized by the BB dielectric function and d = 5.76 Å) levels, as shown in Fig. S14 in the supplementary material. The results reveal a monotonic decrease in both descriptors for all three models, supporting the validity of our full-BEM calculations. In addition, the decay trends for full-ωFQFμ and ωFQFμ-BEM show a one-to-one correspondence, while the full-BEM studies exhibit a smoother decay with distance, further corroborating the robustness of the ωFQFμ-BEM approach.

We conclude this section by remarking on the enhanced computational efficiency of the QM/ωFQFμ-BEM method as compared to the fully atomistic ωFQFμ description (see Table S7 of the supplementary material for the numerical % speed-up). Notably, the computational saving is particularly evident for the largest systems (radius $>$3 nm), for which, as previously discussed, an improved quantitative and qualitative description of the SERS response is also obtained. These findings underscore the potential of ωFQFμ-BEM for its application to large systems, where atomistic details at their surface play a major role in their plasmonic features, offering a substantial reduction in computational cost, up to 58%.

This work has presented a novel multiscale implicit–atomistic approach (ωFQFμ-BEM) to investigate the plasmonic response of metal NPs. ωFQFμ-BEM combines a spherical implicit continuum core embedded within an atomistically defined shell, preserving the essential plasmonic response at the NP surface. The method is theoretically constructed by integrating the fully atomistic ωFQFμ and implicit BEM methodologies for plasmonics and constitutes the first attempt to merge atomistic and continuum descriptions for the study of plasmonic substrates in the context of classical electrodynamics. To validate the ωFQFμ-BEM model, we have initially demonstrated its capacity to reproduce the absorption cross section of spherical Ag nanoparticles, taking the fully atomistic ωFQFμ method as a reference. In particular, we have presented a comprehensive analysis of the different factors influencing ωFQFμ-BEM response: the thickness of the atomistic shell, BEM-ωFQFμ distance, BEM dielectric function, and discretization of the continuum core. Our in-depth analysis suggests that the best results at the implicit continuum–atomistic ωFQFμ-BEM level for Ag NPs are obtained by setting the BEM-ωFQFμ distance to twice the nearest neighbor distance and using the Brendel–Bormann dielectric function while keeping the outer shells described atomistically.

The coupling of ωFQFμ-BEM to a QM region has been presented, giving rise to the QM/ωFQFμ-BEM approach, which is able to describe molecular plasmonics, i.e., how plasmons affect the molecular electronic structure and the resulting spectral signals. To showcase the potentialities of the method, QM/ωFQFμ-BEM has been challenged to reproduce the SERS response of pyridine adsorbed on the tip of Ag Ih NPs. The quality of the method has been assessed by an in-depth comparison with fully atomistic QM/ωFQFμ (reference) and fully implicit QM/BEM spectra. Our results show that, while QM/BEM outcomes strongly differ, QM/ωFQFμ-BEM can reproduce QM/ωFQFμ SERS spectra, specifically for the largest NP sizes. More in detail, an almost perfect agreement of normalized Raman intensities for the largest system is observed, along with a similar qualitative and quantitative outcome for AEF and MEF descriptors with increasing NP radius. As a result of this work, two main points may be highlighted: (i) an implicit description of the NP structure yields substantial discrepancies in simulated SERS spectra as compared to an atomistic description and should be generally avoided; and (ii) it is crucial to describe at a full atomistic level the NP surface, while the core can be safely treated by an implicit description. Finally, the enhanced computational efficiency observed for the largest QM/ωFQFμ-BEM structures yields a computational saving of more than half compared to QM/ωFQFμ. Therefore, large, complex-shaped systems can be afforded. This originates from a substantial reduction of the computational effort required to solve the linear system in Eq. (14), whose dimension is strongly reduced by replacing the core atoms with a continuum description. In addition, memory-related issues are mitigated. To end the presentation, it is worth mentioning that both BEM and ωFQFμ models are formulated in the quasistatic regime, i.e., retardation effects are discarded. Therefore, they appropriately treat systems of dimension significantly smaller than that of the incident wavelength. In addition, in ωFQFμ-BEM, the ωFQFμ and BEM layers interact via electrostatic forces, i.e., charge transfer is inhibited. The above two limitations of the approach will deserve careful investigation in future communications.

See the supplementary material for details on ωFQFμ-BEM equations. Geometrical parameters, PRFs, ρerror, and computational speed-ups. Defined volume for ρerror calculations. Further results on classical and QM/classical responses as a function of core-to-shell distance, dielectric function, and BEM core discretization for spherical and Ih core–shell systems. Graphical representation of pyridine’s normal modes.

The authors acknowledge Professor Stefano Corni (University of Padova) for the fruitful discussions. We gratefully acknowledge the Center for High-Performance Computing (CHPC) at SNS for providing the computational infrastructure. This work has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant Agreement No. 818064).

The authors have no conflicts to disclose.

Pablo Grobas Illobre: Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (equal); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead). Piero Lafiosca: Software (supporting). Luca Bonatti: Methodology (equal); Software (supporting). Tommaso Giovannini: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). Chiara Cappelli: Conceptualization (equal); Funding acquisition (lead); Project administration (lead); Resources (lead); Supervision (lead); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.