Atomistic electrodynamics simulations of bare and ligand-coated nanoparticles in the quantum size regime

The optical properties of metallic nanoparticles with nanometre dimensions exhibit features that cannot be described by classical electrodynamics. In this quantum size regime, the near-field properties are significantly modified and depend strongly on the geometric arrangements. However, simulating realistically sized systems while retaining the atomistic description remains computationally intractable for fully quantum mechanical approaches. Here we introduce an atomistic electrodynamics model where the traditional description of nanoparticles in terms of a macroscopic homogenous dielectric constant is replaced by an atomic representation with dielectric properties that depend on the local chemical environment. This model provides a unified description of bare and ligand-coated nanoparticles, as well as strongly interacting nanoparticle dimer systems. The non-local screening owing to an inhomogeneous ligand layer is shown to drastically modify the near-field properties. This will be important to consider in optimization of plasmonic nanostructures for near-field spectroscopy and sensing applications.


Supplementary Tables
Supplementary Table 1. Morse potential parameters and the scaling factor in electron-density function of EAM potential.

Comparison with GNOR simulations
In Supplementary Fig. 2, we present a comparison between the optical properties of Ag nanoparticle dimers calculated using cd-DIM, DDA and GNOR. The cd-DIM and DDA results are for the nanoparticle dimers discussed in the main text. The GNOR result is based on a 2D model and thus represents infinite Ag nanorods. The GNOR simulations were done in the RF module of the COMSOL Multiphysics package (version 4.4) using the methods developed by Toscano et al 4 . Two 2.46 nm radius rods were placed 1 to 10Å apart (in 1Å step size). The magnitude of the electric field at a point equidistant between the nanorods was calculated as a result from an external field traveling along the long axis of the system, oscillating along the short axis.

Comparison with multilayer Mie simulations
In Supplementary Fig. 3, we present a comparison of the absorption cross section obtained using the multilayer Mie model, cd-DIM and experiment. The Mie results and experimental data are taken from Supplementary Reference 3. Overall, we find good agreement between the different methods. For the larger size we see that the position of resonance peak predicted using Mie theory is in better agreement with the experiments, but the intensity predicted by cd-DIM is in better agreement.

Parametrization of classical force field
The general AMBER force field (GAFF) 5,6 was employed to describe the inter-and intramolecular interactions of the Oleylamine (OAm) ligands. However, there is a lack of the well parameterized force field for describing the Ag-OAm interactions. In order to describe the interactions between the ligand the nanoparticle, we modified a Ag-OAm embedded-atom-method (EAM) potential combined with a Morse potential for describing the Pauli repulsion and DFT-D2 7 dispersion corrections to the long-range van der Waals (vdW) attraction for chemical bonding. This approach has previously been used to simulate ligands binding to Ag surface 8,9 . The total energy of atom i is defined as where F is the embedding energy which is a function of the atomic electron density ρ, and φ is the pair potential interaction between atoms i and j of element types α and β. To describe the chemical bond between Ag-N, we introduce an additional one-way electron-density function to the total electron density 8,9 , which is represented as It is assumed that ρ N−Ag is proportional to ρ Ag−Ag , and, thus we adopt a scaling factor f to evaluate ρ N−Ag : where the pair potential is written in terms of Ag-Ag pair potential taken from the tabulated EAM potential of Ag (http://www.ctcms.nist.gov/potentials) 10 . The Morse potential is given by where D 0,Ag−L , α Ag−L and r 0,Ag−L are the parameters of Morse potential describing the short-range interactions between Ag and OAm ligands. The long-range vdW interaction formulated in the DFT-D2 framework is given by 7 where s 6 is the global scaling parameter taken to be 1.05, and C ij 6 denotes the dispersion coefficient for atom pair ij, which is computed by means of the combination rule, In Eq. 5, f damp (r ij ) denotes the Fermi-type damping function.
Typically d=20 and R ij 0 denotes the sum of the vdW radii of atom pair ij. The reference data used for the parametrization of the EAM-like potential were generated from DFT calculations, where we chose a model system consisting of n-butylamine (BAm) and a Ag 56 cluster, see Supplementary Figure 6. We considered the perpendicular and parallel orientations of the BAm molecule with respect to the Ag(111) facet. For the DFT calculations we used the Becke-Perdew (BP86) XC-potential 12,13 with a triple-ζ polarized Slater type (TZP) basis sets to optimize the geometries using the ADF program package (http://www.scm.com) 14 . The 1s-3d core was frozen for Ag and scalar relativistic effects were taken into account by using the zeroth-order regular approximation (ZORA) [15][16][17] . Moreover, the dispersion energy was calculated by the DFT-D3 approach 18 . Only small differences in the geometry optimization were found by the DFT-D2 and DFT-D3 and thus the simpler D2 dispersion model was used in building the force field.
The parameters in the Morse potential and the scaling factor in the one-way electron-density function were found by minimizing the following cost function Here f ijk represents the force along the Cartesian component k originating from the Ag cluster on atom j in the molecular configuration i. E inter i denotes the interaction energy between BAm and Ag cluster in configuration i. The optimized parameters are collected in Supplementary Table 1. 16 grid points used in the parametrization were generated by perpendicularly translating the BAm molecule around the optimized geometry in step of 0.1Å with respective to the Ag(111) facet, while the geometries of BAm and Ag cluster are fixed.
A comparison of interaction energies obtained from DFT and the force field for BAm is represented in Supplementary Fig. 5 and the structure difference is illustrated in Supplementary Fig. 6. Overall, we see that the binding energy for the two orientations is correctly described using the force field. Comparing the geometries optimized by the force field and DFT for the perpendicular and parallel configurations, we find the RMSD are 0.042Å and 0.439Å, respectively.

Parameterization of atomic radii to describe the static polarizability
The atomic radii representing the surface and bulk atoms in cd-DIM were chosen to describe the mean static polarizability of small metal clusters as compared with DFT results. The size-dependence of the static polarizability can be described by a jellium model 1,19 .
where r WS is the Wigner-Seitz radius of the bulk metal and δ represents the spillout of the electrons from the surface of a metallic sphere. In Supplementary Fig. 1 we present the data obtained using cd-DIM for Ag clusters with N ≤ 68 that was previously used to study the polarizability using the capacitance-polarizability interaction model (CPIM) 1  To obtain the atomic radii needed for describing the ligands we choose to minimize the difference between cd-DIM and DFT polarizability. The multiple configurations of BAm were selected from MD trajectories and used to calculate the polarizability using DFT. In the MD simulations, BAm was placed in the center of 15×15×15 A 3 periodic box, and the linear and angular momenta were zeroed. Followed by energy minimization, the temperature was increased gradually from 0 to 300 K under the NVT ensemble within 50 ps, and subsequently the whole system was equilibrated for 10 ns with an integration time step of 0.5 fs. Periodic boundary conditions with a cutoff distance of 12Å for vdW and electrostatic interactions were employed. The static polarizability was calculated by DFT at the B3LYP/aug-cc-pVDZ level of theory using the program package NWChem 6.1.1 22 .
The atomic radii of BAm were then determined by minimize the following cost function where α i,αβ represents the each component in the polarizability tensor of configuration i, and N is the number of configurations, which is set to 50. The comparison between the polarizability calculated by DIM and DFT is shown in Supplementary Fig. 4.