Abstract

Molecular dynamics (MD) is a computer simulation technique that studies particle models by stimulating the time evolution of a system of interacting particles. MD is commonly employed in materials science, nanotechnology and biochemistry to study various processes at the atomic scale. The techniques of MD have also been adapted for use in fluid dynamics and astrophysics. In this article, we introduce the basic science and computational methods behind MD and explore some of its current applications.

Given a system consisting of n particles subjected to interactions described by a given force field, can the trajectory of every particle be predicted? 18^th century French astronomer and mathematician Pierre-Simon Laplace said of solving the above problem, called the n-body problem, “Given for one instant an intelligence which could comprehend all the force by which nature is animated and the respective situation of the beings who compose it—an intelligence sufficiently vast to submit these data to analysis—it would embrace in the same formula the movements of the greatest bodies of the universe and those of the lightest atom; for it, nothing would be uncertain and the future, as the past, would be present to its eyes.”

Molecular dynamics (MD) solves the n-body problem by numerically integrating the equations of motion for every particle. There are two common forms of MD, ab-initio MD and classical MD. Ab-initio MD solves the many-body problem for a system governed by the Schrödinger equation with the Born-Oppenheimer approximation (without this approximation, the Schrödinger equation cannot be solved even numerically for many systems). At every timestep, the Schrödinger equation for electrons is solved to compute the potential. The position and velocity of every nucleus are then propagated by numerically integrating the classical equations of motion. To solve the Schrödinger equation for electrons, various approximations derived with methods such as the tight-binding model or density functional theory are used. Ab-initio MD can simulate quantum mechanical phenomena, such as bonding between atoms and heat conduction in metals, but the computational cost limits the size of systems that can be studied with this model to a few thousand atoms.

Classical MD simplifies the ab-initio approach by replacing the quantum mechanical potential computations with a parameterized empirical potential function that depends only on the position of nuclei. This potential function is either fitted to experimental data or derived from the underlying quantum mechanical principles. Classical MD can be used to study systems with up to millions of atoms but cannot simulate phenomena that depend on electron behavior. We will restrict our discussion to classical MD.

The name molecular dynamics is misleading because given the almost universal nature of the many-body problem, MD can model any system that involves interaction between particles. MD is employed in materials science to study a variety of phenomena, including the behavior of materials under stress, crack propagation, and how various defects affect the strength of materials. MD is also commonly used in biochemistry to study the behavior of macromolecules. The techniques of MD have also been generalized to thermal science and astrophysics, where particle models are used to study various hydrodynamic instabilities and phase transitions in thermal science and the structure of the universe in astrophysics.²

Interaction computations

The computation of interactions between atoms is typically the most computationally expensive segment of a MD simulation. For a pair-wise potential (a potential described by a potential function that does not account for the environment of atoms), the interaction computations scale with order O(N²). Various methods have been developed that reduce the order to O(N) for short-range potentials (potentials that decay faster with respect to distance than the dimension of the system are considered short-range) and O(NlogN) for long-range potentials. The most common methods are the linked-cell and neighbor list methods for short range potentials and Ewald summation method for long-range potentials.

Both the linked-cell and neighbor-list method assume that the potential and its spatial derivative at the location of an atom can be approximated by the superposition of the potential contributions from atoms within a certain interaction cutoff distance. Mathematically, the potential function U(r_ij) can be approximated by the truncated potential function Û(r_ij).

But to compute the distance between every atom pair costs O(N²) calculations. The linked-list method decomposes the system spatially into cells of side length greater than or equal to the cutoff distance. Every atom only interacts with atoms in the same cell or an adjacent cell so only distances between atom pairs in the same and adjacent cells must be computed. Atoms are resorted into cells every 10 to 20 timesteps. This reduces the complexity of the computation to O(N) .²

The linked-list method is commonly used in parallel molecular dynamics codes because it simplifies communications between threads. Every thread is assigned a cell and is passed a list of the atoms in its cell and every adjacent cell every time atoms are resorted into cells.

Figure 1 Linked-cell method: To compute the interactions of the filled in atom, only atoms in the same (dark-shaded) and adjacent (light-shaded) cells are considered. Only atoms within the shaded disk interact with the darkened atom.

The neighbor list method creates lists of all atoms within a certain distance of every atom so every atom only interacts with atoms in its neighbor list. The neighbor list cutoff distance must be greater than the interaction cutoff distance and the difference determines the frequency at which the neighbor lists must be recreated. The neighbor list method also reduces the complexity of the computation to O(N) (2).

Figure 2 Neighbor list method: The neighbor list for the filled-in atom contains all the atoms in the light-shaded disk. Only atoms within the dark-shaded disk interact with the filled in atom.

For long-range potentials, the interaction cutoff distance must be very large for the potential at the location of an atom to be approximated by the superposition of the potential contributions from atoms within the cutoff distance. The Ewald summation method is the most common method and reduces the computational complexity to O(Nlog(N)) ³ by decomposing the system spatially into cells and assuming the interactions of an atom with every other atom in the system can be approximated by its interactions with every other atom in the same cell and their images in periodic cells tessellated along all dimensions.

Mathematically, this can be expressed as:

where N is the number of atoms in the same cell, n_x, n_y, n_z are half the number of periodic cells in the x, y, z directions respectively, Ø(r_ij) is the part of the function that does not decay with respect to distance, L is the side length of the cell and d is the dimension of the system. The Ewald summation method then decomposes the above sum into a short-range term that can be summed using techniques for short-range potentials and a long-range term that can be summed with Fourier methods. A detailed mathematical introduction to the Ewald summation method is beyond the scope of this article and we refer the reader to ³ for a detailed introduction to the Ewald summation method.

Integration methods

The solution to an n-body problem can be chaotic and it is fruitless to seek an accurate trajectory beyond several “collision steps”. The main criteria for choosing an integration method for the equations of motion in MD are therefore energy conservation, the ability to accurately reproduce thermodynamic ensemble properties and computational requirement. Interaction computations are the most computationally expensive segment of an MD simulation and an integration method requiring more than 1 interaction computation per timestep is inefficient unless it can maintain the same accuracy with a proportionally larger timestep. This criterion disqualifies the common Runge-Kutta methods.

Two commonly used families of integration methods are the Verlet algorithm and its variations and predictor-corrector methods. The Verlet algorithm can be derived from the Taylor expansion of the position variable (r):

Adding the expressions for r(t+h) and r(t-h) and solving for r(t+h):

One drawback of the Verlet algorithm is the lack of a velocity term in the expressions, making it difficult to compute the atomic velocities. Common variations of the Verlet algorithm designed to correct this drawback are the leap-frog algorithm and velocity-Verlet algorithm.

The predictor-corrector family of methods (PC methods) all consists of 3 steps. First, atomic positions and velocities are propagated according to the Taylor expansions of the position and velocity variables. Then, the interactions are computed based on the new positions and the accelerations are compared to the accelerations predicted with the Taylor expansion of the acceleration variable. Finally, the positions and velocities are “corrected” according to the difference between the computed and predicted accelerations. Although PC methods require 2 interaction computations per timestep, they can maintain the same accuracy with a timestep that is more than twice as long.³

PC methods have mostly been discarded in favor of the simpler Verlet methods because Verlet methods give better energy conservation and are easier to implement. We refer the reader interested in a derivation of the basic PC algorithm to. ¹

Constraint dynamics

Geometric constraints can be implemented in the context of the above integration methods to simulate rigid molecules. The rigid molecule approach is not appropriate for large molecules. Methods to simulate flexible molecules have been developed although we will not introduce them in this article. We refer the reader to¹ for a detailed introduction to the simulation of flexible molecules. The most common method is the SHAKE algorithm and its variations. The SHAKE algorithm first advances the atomic positions and velocities ignoring the constraints. The new positions and velocities are then corrected by introducing a constraint force, the magnitude of which is determined by the displacement of the new positions from their constrained positions.

The k^th constraint constraining the distance between two atoms can be expressed mathematically as:

where r_i and r_j are the positions of atoms i and j and d_k is the bond length between them. These constraint equations represent forces that are described by their action on the dynamics of the system and are not derived from the potential. These forces can be expressed within the equations of motion as constraint potentials:

where n is the number of constraints and λ is a constant that must be computed in order for the constraint to be satisfied. We numerically integrate the equations of motion with the Verlet algorithm.

where is the unconstrained position of atom i. r_i and r_j must obey the constraint equations σ_k:

where r_k is the constrained bond length. The value of λ_k can be solved for with standard root-finding techniques. In the SHAKE algorithm, the system of k equations is solved with iterative methods [1].

Thermostats and barostats

MD with the classical equations of motion simulates a system at constant volume and energy. Every timestep, a new state of the system, called a microstate, is generated. Because every microstate generated has the same volume and energy, MD with the classical equations of motion can be thought of as sampling from a set of microstates with a given number of atoms, volume and energy, called a microcanonical/NVE ensemble (so called because the Number of particles, Volume and Energy are constant).

It is often desirable to sample from a canonical/NVT (Number of particles, Volume and Temperature are constant) or isothermal-isobaric/NPT (Number of particles, Pressure and Temperature are constant) ensemble to maintain compatibility with experimental studies. Various methods have been developed to constrain temperature and pressure and are called thermostats and barostats respectively. Most modify the system dynamics by modifying the equations of motion. Common thermostat algorithms include the Berendsen thermostat and Nosé-Hoover thermostat and common barostats include the Parrinello-Rahman method.

The equations of motion for most thermostat algorithms take the general form:

where γ can be thought of as a coefficient that adjusts the system kinetic energy to match the desired temperature. The effect of a thermostat algorithm on system dynamics can also be described as a scaling of internal velocities (velocities relative to center of mass) at every timestep by a factor λ. Mathematically, this can be expressed as:

where is the unscaled internal velocity and h is the timestep.

Because the temperature cannot change instantaneously, λ(t,0) = 1.

If we let in (*), then the above expression agrees with (*).

The Berendsen thermostat is based on Newton’s law of cooling:

where T₀ is the heat bath temperature and τ is an empirical coupling parameter that quantifies the strength of the coupling between the system and heat bath. The finite difference quotient is:

The velocity scaling factor at every timestep is determined by the ratio between the system temperature and heat bath temperature and the coupling parameter.

Note that the coupling parameter affects the system dynamics like a damping parameter. The Berendsen thermostat, therefore, does not generate a true canonical velocity distribution because the coupling parameter damps out system temperature fluctuations.

Both the Nosé-Hoover themostat and the Parrinello-Rahman method are based on the extended system approach. In this approach, additional virtual degrees of freedom are introduced into the system that can be used to regulate system properties.

A detailed introduction to the Nosé-Hoover themostat and the Parrinello-Rahman method is beyond the scope of this article. We refer the reader ⁴ for a detailed introduction to the Nosé-Hoover thermostat and²
for a detailed introduction to the Parrinello-Rahman method.

Conclusion

According to American physicist Richard Feynman, “all things are made of atoms and everything that living things do can be understood in terms of the jiggling and wiggling of atoms”. MD is a computer simulation technique that simulates the “jiggling and wiggling of atoms” and has been hailed as a way of predicting the future by animating nature’s forces. In this article, we introduced the basic method behind MD to provide the reader with the necessary background knowledge to understand an MD simulation.

References

[1]. Rapaport, D. C.The Art of Molecular Dynamics Simulations. Cambridge : Cambridge University Press, 2004.
[2]. Griebel, Michael, Knapek, Stephan and Zumbusch, Gerhard. Numerical Simulation in Molecular Dynamics: Numerics, Algorithms, Parallelization, Applications. Berlin : Springer, 2007.
[3]. Leach, A. R.Molecular Modeling: Principles and Applications. Upper Saddle River : Prentice Hall, 2001.
[4]. Frenkel, Daan and Smit, Berend.Understanding Molecular Simulation: From Algorithms to Applications. San Diego : Academic Press, 2002.
[5]. Hünenberger, Philippe. Thermostat Algorithms for Molecular Dynamics Simulations. Advanced Computer Simulation. Berlin : Springer, 2005, pp. 105-149.

Abstract

One of the most significant developments in the field of Raman spectroscopy was the discovery of Surface-Enhanced Raman scattering (SERS). Raman signals from the early days of SERS were known to be enhanced by a factor up to a million-fold, and more recently by extraordinary factors of up to 10^12-14. To understand what this means we should understand how Raman scattering works, and then apply this understanding to the Surface-Enhanced case. SERS can be thought of as at least two kinds of enhancement: an electromagnetic enhancement which is well understood and a chemical enhancement which is not completely understood. Charge-Transfer has been shown to contribute some chemical enhancement, but more mysteries remain. We are working to characterize the wavelength dependence of the Anti-Stokes to Stokes ratio, thus deriving what contributions are made by metal-adsorbate-complex transitions approximately resonant with both the laser and surface plasmon frequencies in the absence of Charge-Transfer effects. Work by this group and by others indicates that anti-Stokes to Stokes ratios can be sensitive indicators of this chemical enhancement mechanism in SERS when it occurs.

Background: Raman Scattering

A new type of light-matter interaction was experimentally verified by C.V. Raman and K.S. Krishnan in 1928, now called Raman scattering.¹ In contrast with Rayleigh scattering, it is the scattering of photons with a change in frequency. Though the effect holds for all wavelengths of light, it is easier to consider a monochromatic case. During a Raman event, an incident photon of a particular wavelength and energy interacts with the vibronic modes of a matter system. The term vibronic indicates a combination of electronic states and vibrational states. The energy of the system increases by an amount allowed by vibronic transitions and a photon having lost this exact energy is re-emitted. The molecule thus moves from a lower energy state to a higher energy state, while photons of a higher energy are absorbed and photons of a lower energy are emitted. This is called a Stokes shift, but each vibronic transition also has an Anti-Stokes shift. (See Figure 1.) In the Anti-Stokes case, the matter is already in an excited vibrational state. An incident photon interacts with the excited system, resulting in the emission of a higher energy photon and a loss of energy in the system as it returns to a lower energetic state. This happens less often and is dependent on how many molecules in the system are in an excited state, generally due to thermal interactions. The probability of any Raman event occurring is on the order of one for every million incident photons, but it is highly dependent on the polarizability of the matter, and the frequency of the incident radiation. The ratio of Anti-Stokes events to Stokes events is a useful quantity for observing resonance between vibronic transitions and incident radiation.

Since the advent of laser technology in the mid-20th-century, scientists have had access to collimated beams of nearly monochromatic radiation. From this, Raman spectroscopy was developed to probe the vibronic transitions of matter systems by examining the emitted Raman light with a spectrometer. If the energy of the incident radiation is known, the energy of the emitted radiation can be measured, and the energy gained by the system can be calculated exactly. This energy then corresponds to vibronic modes in the system, offering experimental verification of theoretical calculations. Like Infrared or Fluorescent Spectroscopy, Raman Spectroscopy identifies groups of molecules by their chemical bonds, but also provides additional information. The vibrational modes of matter are easily characterized and studied, and a large body of work dedicated to this field exists. In Raman spectroscopy, vibrational modes are called Raman modes. To display Raman spectra, the wavenumber of the incident radiation (usually a laser) is marked as zero. In what is known as the wavelength shift, the wavenumber of each Raman mode is shifted from zero, corresponding to the energy of that mode. For example, a carbon-carbon stretch occurs at approximately 1600 wavenumber shift. It will always occur at 1600 wavenumber shift, regardless of the frequency of the incident radiation. Changing the incident radiation changes the absolute wavenumber position, but not the relative shift. Note that varying the surrounding environment or the physical state of the matter may shift the Raman modes slightly. Figure 4 shows two Surface-Enhanced Raman spectra of 1-dodecanethiol. Several carbon-carbon modes are visible between 1050 and 1200 wavenumber shift.

Background: Surface-Enhanced Raman Scattering

Due to the rarity of Raman events, the Raman signal is often very weak relative to the intensity of incident radiation. The experimentalist desires very intense incident radiation to be ensured of a strong signal. However, over-heating the sample can lead to unwanted side effects, so a balance must be sought between strength and authenticity of the signal. One method to mitigate this is to place the matter in close proximity with a noble metal, such as gold, silver, or copper. Because the intensity of Raman events is dependent on the polarizability of the matter, the strength of the signal is closely dependent on the number of free electrons in the system. Noble metals have an abundance of free electrons that become polarized in resonance with the incident radiation and the molecules being studied, resulting in an increased intensity of Raman photons. This is called Surface-Enhanced Raman Scattering, or SERS. When light interacts with a noble metal, the free electrons on the surface tend to polarize in resonance with the light, and this collective motion is called a surface plasmon. The peak resonance of the surface plasmon depends on the type of noble metal and its geometry. Understanding surface plasmons and their wavelength dependence is critical to understand surface enhancement.

In the late seventies, electrochemists took spectra of pyridine adsorbed on silver, copper or gold electrodes electrolyte solutions. ^3-5 Scanning the potential of the electrode (by comparison with Saturated Calomel Electrode), they demonstrated variations in the Raman signal. It was discovered that oxidation/reduction cycling would dissolve and redeposit monolayers of metal, augmenting the SERS enhancement by two orders of magnitude. This additional enhancement was attributed to molecular roughening of the surface of the electrodes.⁶ The SERS intensity was found to increase with the number of deposited monolayers for small numbers of layers.⁷ Allen et al. were able to show that multiple mechanisms contribute specified magnitudes to overall enhancement. They suggested that approximately 3.5 x 10¹ enhancement could be contributed from a surface-roughness mechanism, and 2 x 10³ enhancement could come from roughness-independent mechanisms for pyridine adsorbed on copper electrodes.⁸

Wavelength Dependence of SERS: Multiple Mechanisms

Pockrand et al. determined that the multiple mechanisms were not limited to substrate activity. They collected Raman spectra of pyridine on evaporated silver films and on silver electrodes, and compared SERS-active films with SERS-inactive films investigating wavelength dependence of Raman intensities. ⁹ They explicitly noted that excitation profiles could be plotted showing SERS intensities dependent on the laser excitation wavelength, with a single maximum in the visible region for each vibrational mode. (See Figure 2.) Their paper demonstrated that all the excitation profiles they took have similar resonance behavior, but that the curves shift to greater energy of the incident light with the increasing energy of the vibration. The surface plasmon enhancement of silver was qualitatively the same with evaporated silver films and with silver electrodes, ruling out mechanisms solely dependent on the substrate. Regarding the adsorbed molecule, they found that doublet signals were obtained from thick layers of pyridine on SERS-active films; one of the pair correlated to the bulk pyridine and the other to a red-shifted version of the same vibrational mode in the surface pyridine. This corroborates the conclusion that the SERS response for pyridine is dependent on both surface plasmon activity attributed to the silver substrate and excitation energy of the laser. Furthermore, the large enhancement factors from SERS-active films could not be described by local field effects alone, and Pockrand concluded that a chemical contribution must play a role.

The multiple mechanism picture has prevailed over the years, leading to the proposal of various models, more or less dependent on specific data obtained from SERS experiments. In general, it is now well-accepted that a local-field electromagnetic enhancement occurs, which resonates with vibronic modes in the molecule and thus substantially increases the intensity of Raman scattered radiation. This mechanism requires that the molecule’s Raman-shifted vibronic frequencies fall within the spectral width of the substrate surface plasmons, although direct proximity between metal and molecule are not required. ^{10, 11} It is also well-accepted that this mechanism does not completely describe SERS. For physisorbs/chemisorbs, molecule-specific chemical effects are found in enhancement patterns and intensities. ¹⁰ A Charge-Transfer mechanism, requiring direct proximity between metal and molecule, has been suggested to explain these.^11-13 Charge-Transfer involves the overlapping of molecular orbitals to provide pathways for the migration of electrons from the metal to the adsorbate, or vice versa. This is thought to provide enhancement by modifying the Raman polarizability tensor, causing increased intensity in Raman scattered radiation, though some scientists argue that the Charge-Transfer mechanism should be regarded as an altered adsorbate complex rather than an enhancement mechanism.¹¹

Early work with electrodes, colloids and island films was insufficient to fully understand the complexity of SERS substrate plasmon resonance activity. Researchers called for experiments with increased control of surface roughness, but technology was lacking. ⁸ A solution to this problem has been developed at Rice University and involves the fabrication of nanoscale particles, like nanoshells. The study of nanoparticles has led to major advancements in the understanding of SERS, such as Jackson et al. recently investigating the relationship between nanoshell geometry and SERS enhancement. ¹⁴ Using nanoshells, a new monodispersed substrate, Jackson et al. showed that, in particular for those modes which the laser is resonant with, the Raman Effect is enhanced proportionately with the density of the nanoshells. ¹⁴ In other work by the same authors, it was discovered that SERS intensity was linearly dependent on the density of the nanoshell substrate, proving that the SERS response in their experiments was due to single nanoshell resonance response, not dimer resonances. ¹⁵ It is important to note that tuning the laser wavelength to excite single nanoshell resonances or dimer resonances results in data corroborating one theory or the other. This reinforces the notion that wavelength dependence is critical to understanding SERS mechanisms.

Maher et al. specifically investigated wavelength dependence of SERS on multiple substrates by plotting the anti-Stokes to Stokes ratio for p-mercaptoaniline (pMA) at different excitation energies and discovered that asymmetries arise. ¹⁶ At a fixed temperature a calculable amount of anti-Stokes emission is usually expected based on the thermal activity of the sample. They used the ratio of anti-Stokes to Stokes intensities as a sensitive indicator of resonant conditions involving both the metal and the adsorbate molecule. They argued that this asymmetry is systematically structured according to underlying resonances that need more than a single wavelength to be seen. One series of data they took for laser line 633 nm produced a higher anti-Stokes to Stokes ratio than expected, and another series using a 514 nm excitation produced a lower ratio than expected. (See Figure 3.) This may correspond to an anti-Stokes emission intensity peak around the 585 nm to 610 nm region. Their data lacked the excitation energy density necessary to determine this with certainty, but they did show that the anti-Stokes to Stokes ratio is a potentially valuable means of mapping the resonant interaction between the probes and the metal.

Advancing Our Understanding

We anticipate that in some cases unexpected resonances (called “hidden” by Maher et al. ¹⁶) for the metal-molecule-complex will be observed, even in the absence of Charge-Transfer processes. We are specifically interested in 1-dodecanethiol, though we also intend to do experiments with p-mercaptoaniline, due to the other work that has been conducted at Rice University. ^10,16 Initial observations for the 1-dodecanethiol and gold colloids system using a Renishaw Invia Raman Microscope have been made. (See Figure 4.) We are currently building a Tunable Raman Darkfield Microscope system to restrict the region of measurements to those where the gold colloids have deposited in a smooth and dense manner. Using a dye laser, we can tune continuously from 558 nm to 606 nm, providing good coverage of the visible spectrum in the region where gold colloids enhance strongly. Compiling excitation profiles as a function of excitation energy, we will calculate the theoretical thermally-allowed anti-Stokes to Stokes ratio for the sample and then examine deviations from this result. Building on the work of Gibson et al. (see subsequent section of this paper), ¹⁰ we may prove or disprove the theory of hidden resonances for the metal-molecule-complex in the absence of Charge-Transfer for pMA attached to gold nanospheres and nanoshells. This work can later be applied to other substrates and molecules to provide data to support the general case.

Theoretical Support

As theoretical support for our work, we turn to Gibson et al. who independently ran theoretical density-matrix calculations for pMA and obtained results corroborating the notion that unexpected resonant behaviors between molecule-metal-complex may play a significant role in Surface Enhancement. ₁₀ Their model was able to obtain anti-Stokes resonant excitation behavior in the absence of Charge-Transfer, supporting our hypothesis that these unexpected resonances occur when the anti-Stokes emission is in resonance with vibronic transitions of the molecule-metal-complex. Gibson et al. suggest that anti-Stokes emissions which are resonant with a metal-adsorbate molecule may provide anti-Stokes intensity peaks which then fall away in intensity on either side of the peak. Their calculations agree with the work of Maher ¹⁶ and may explain why his observed anti-Stokes to Stokes ratio is asymmetrical to thermal expectations; lower or higher, depending on laser excitation energy.

Acknowledgements

The author would like to thank Rice University and the Rice Quantum Institute. He also acknowledges funding provided by the National Science Foundation and hopes that Research Experience for Undergraduates programs will continue to be funded for the long-term development of research sciences.

References

1. Long, D.A. “Raman Spectroscopy”, (McGraw-Hill, Great Britain, 1977).
2. Image from Vrije University Brussels, 2007, Dec 18, . Accessed 2008 Mar10. Axis labels added.
3. M. Fleischmann, P.J. Hendra and A.J. McQuillan, Chem. Phys. Letters 26, (1974) 123.
4. U. Wenning, B. Pettinger, H. Wetzel, Chem. Phys. Lett. 70, (1980) 49.
5. J. Billmann, G. Kovacs and A. Otto, Surface Sci. 92, (1980) 153.
6. D.L. Jeanmarie and R.P. Van Duyne, J. Electroanal. Chem. 84, (1977) 1.
7. B. Pettinger and U. Wenning, Chem. Phys. Letters 56, (1978) 253.
8. Craig S. Allen, George C. Schatz, Richard P. Van Duyne, Chem. Phys. Letters 75, (1980) 201.
9. I. Pockrand, J. Billman, A. Otto, J. Chem. Phys. 78 (11): 6384-6390 1983.
10. J.W. Gibson, B.R. Johnson, J. Chem. Phys. 124, 064701 (2006).
11. M. Moskovits, Rev. Mod. Phys. 57 (3): 783 1985.
12. M. Osawa, N. Matsuda, K. Yoshii, I. Uchida, J. Phys. Chem. 1994, 98 12702-12707.
13. J.R. Lombardi, R.L. Birke, T. Lu, J. Xu, J. Chem. Phys. 84 (8) 4174, 1986
14. J.B. Jackson, S.L. Westcott, L.R. Hirsch, J.L. West, N.J. Halas, Appl. Phys. Lett. 82, 257 (2002).
15. J.B. Jackson, N.J. Halas, Proc. Nat. Acad. Sci. USA 101, 17 930 (2004).
16. R.C. Maher, J. Hou, L.F. Cohen, E.C. Le Ru, J.M. Hadfield, J.E. Harvey, P.G. Etchegoin, F.M. Liu, M. Green, R.J.C. Brown, M.J.T. Milton, J. Chem. Phys. 123, 084702 (2005).