For those of us working with nanoparticles, it became a common task to simulate the plasmon peak position as a function of geometry. Even if most papers claim it is an easy task, we all know that it is easy once you know how to do it; this article therefore aims to help those who are starting with simulations using the excellent ADDA (Amsterdam-Discrete Dipole Approximation) package, aiming to have spectra of gold nanoparticles. I will focus on gold nanorods, but the examples provided below should be easily extendable to other geometries.
Installing the ADDA package has no misteries; it is precompiled for Windows and you can download the source code for linux (in which I’m going to base the following examples.) As far as I understand, there is no easy way of generating an absorption spectra but to calculate the cross sections at different wavelengths. This of course, is a tedious task that can be easily performed by scripting via Python (of course almost any other languages, including Matlab, may work). I will stick to Python 3+, with basic packages as numpy and scipy for some math operations. First things first, let’s install and test ADDA (of course, if you are on Linux, Windows users may skip this paragraph). Download and extract the zip file. Go into the src folder and run
if you want the sequential execution of the program (if you have a multi core processor this is not the best way). You can also run
to compile the parallel version of it (using all the cores). There is also an option for running on GPU but I couldn’t try it. Once compiled, you will see the presence of executables called adda_seq and/or adda_mpi inside their respective folders. I recommend linking them in your /usr/bin/ folder, in order to run the with a simple command. To test the installation, just run adda_seq (or adda_mpi) and check the output (it just computes the cross sections of a gold nano-sphere). Now that everything is running, let’s start with the interesting part. The first step is to get the values for the gold dielectric function. You can check the papers you want and extract the values from them, in my case I use the Christy and Johnson (a reference in the field). Once you have the values, you need a function to interpolate them; I’m using cubic splines (check that your simulation wavelength is not falling outside the given values and the objections that appeared to those reported values.) The code for that is quite plain and you can see it below. Check the names I’ve given to the file with the dielectric values. You can fetch all the codes explained here from my GitHub repository
import numpy as np from scipy import interpolate # For spline interpolation class Au_index(): def __init__(self): #Read the file data = np.loadtxt('METALS_Gold_Johnson.txt',skiprows=1,unpack=True) wavelength = data real_index = data imag_index = data # Spline interpolation of real and imaginary parts ind = wavelength.argsort() wavelength = wavelength[ind] real_index = real_index[ind] imag_index = imag_index[ind] self.real_tck = interpolate.splrep(wavelength,real_index) self.imag_tck = interpolate.splrep(wavelength,imag_index) def r_i(self,wavelength): """ Returns the real and imaginary values of the refractive index evaluated in the desired wavelength. """ real = float(interpolate.splev(wavelength,self.real_tck)) imaginary = float(interpolate.splev(wavelength,self.imag_tck)) return complex(real,imaginary) class H2O_index(): def __init__(self): #Read the file data = np.loadtxt('LIQUIDS_Water_Hale.txt',skiprows=1,unpack=True) wavelength = data real_index = data imag_index = data # Spline interpolation of real and imaginary parts ind = wavelength.argsort() wavelength = wavelength[ind] real_index = real_index[ind] imag_index = imag_index[ind] self.real_tck = interpolate.splrep(wavelength,real_index) self.imag_tck = interpolate.splrep(wavelength,imag_index) def r_i(self,wavelength): """ Returns the real and imaginary values of the refractive index evaluated in the desired wavelength. """ real = float(interpolate.splev(wavelength,self.real_tck)) imaginary = float(interpolate.splev(wavelength,self.imag_tck)) return complex(real,imaginary)
You can see that I’ve defined two classes, one for gold and one for water; they both work in the same way, returning a complex value for the dielectric constant at a given wavelength. IN the future, the water class will not be used, but you can keep it for reference as well. The class doesn’t do any error checking (for instance if you set a wavelength outside the range) and therefore you have to be fully aware of what you are doing. For the time being we are going to assume that water has a constant dielectric value of 1.33. And yes, there is no need of defining two different classes, one could have just defined a ‘material’ class and just change the file being read. The way ADDA works needs to normalize everything by the dielectric constant of the medium, in this way the wavelength of the incident beam is going to be divided by 1.33, etc. For simple geometries (like a sphere or a rod) ADDA comes with a really handy function to determine geometries. You should check the documentation for further information. In my case, I represent gold nanorods as capsules (a cylinder capped by two hemi-spheres.) The parameters passed to the cylinder are the aspect ratio (careful with the definition, since it is not length per radius!) and the dimension (the length). For rods the orientation is important (since the incident beam will be polarized, you want to be sure the rod is parallel to the polarization.) So far, for simulating a 25×60 rod and wavelength of 650nm you can run the following code:
adda -shape capsule 1.4 -size 0.025 -grid 16 -orient 0 90 0 -lambda 0.489 -m 0.1136 2.707
Remember that the way of executing adda will depend on how you installed it in your computer and it is your task to find it out. The output should be something similar to this (the last few lines):
Cext = 0.02126180812 Qext = 20.37303979 Cabs = 0.01680466156 Qabs = 16.10220714
And voilà, you have calculated the cross-section of a rod for that particular wavelength. Check the last argument in the simulation. It is the dielectric constant passed as two numbers. Of course we can add the class we defined earlier for calculating the values at different wavelengths. We can also include everything in a for loop to generate the full spectrum in a range of wavelengths, we’ll see how to achieve this in the following pages.