# Fitting Complex Metal Dielectric Functions with Differential Evolution Method

Posted on Wed 13 April 2016 in Plasmonics

The real and imaginary part of dielectric permittivity of the metals is important to simulate the optical properties of metal films and nanoparticles. Permittivity data is obtained experimentally by ellipsometry and is fitted with analytical models. The most common model for fitting experimental data is with Drude-Lorentz model shown below. \$\$\\epsilon(\\omega)=1-\\frac{f\_1\\omega\_p\^2}{(\\omega\^2+i\\Gamma\_1\\omega)}+\\sum\_{j=2}\^{n}\\frac{f\_j\\omega\_p\^2}{(\\omega\_{o,j}\^2-\\omega\^2-i\\Gamma\_j\\omega)}\$\$ The first term is the Drude part. It represents the response of electron in the Fermi sea/conduction band when it sees external oscillating electric field (these transitions are called as intraband transitions). The Drude term has the plasma frequency (\$\\omega\_p\$, oscillator strengh (\$f\_1\$) and damping term (\$\\Gamma\_1\$). The rest of the terms represent the Lorentz oscillators with specific resonance frequencies (\$\\omega\_{o,j}\$), oscillator strengths (\$f\_{j}\$) and damping terms(\$\\Gamma\_j\$) associated with them. They represent electron excitation from one band to another following an external oscillating electric field at certain resonance frequencies (these transitions are called as interband transitions). Just out of curiosity, I wanted to do this fitting by myself for a long time now. So I attempted to fit the data with the commonly used fitting method, the least squares method, which involves minimization of the square of the error between the model and experimental data, while changing the parameters. However I was not getting good fitting results even with several attempts, this may be due to the following problems with this method: - The fits were very sensitive to the initial guess values of the parameters (oscillator strengths, damping terms, etc). - The method is a local minimization technique. What that means is that the least square algorithm (Levenberg - Marquardt ) is changing the parameters in parameter space and searching for the minimum of sum of squares of difference between model and data. While doing this it can get trapped in a local minimum and report the parameters at that point as if it was as global minimum. For more information on the problems of least-square method. See the last section of the article [here](http://www.itl.nist.gov/div898/handbook/pmd/section1/pmd141.htm). To avoid these problems, researchers use global optimization techniques. See, [Rakic et al.](https://www.osapublishing.org/ao/abstract.cfm?uri=ao-37-22-5271%20) and [Djurisic et al.](%20http://iopscience.iop.org/article/10.1088/1464-4258/2/5/318/meta). Rakic et al. used simulated annealing method , whic is a global optimization method and does a great job fitting experimental data. I wanted to use [Scipy's simulated annealing method](http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.optimize.anneal.html) to fit the dielectric functions, but found that simulated annealing algorithm is deprecated in scipy 0.14 and they suggest using basin-hopping or differential evolution (DE) algorithms instead. I gave DE method a try because I didnot see any research article using DE algorithm for fitting dielectric data. If any of the readers find an article that uses DE to fit complex dielectric function, please point me to the article. DE algorithm is a stochastic method that does not involve finding a gradient but rather involves mutation and recombination of random population of solutions. I dont completely understand the guts of the algorithm, but you can find more information [here](https://en.wikipedia.org/wiki/Differential_evolution). Below is my code of how I do the fitting of gold dielectric function with python and lmfit module. lmfit is a python module that provides a better and convinient fitting and optimization interface for Scipy optimize/minimze methods. The default minimization technique for lmfit is least square (Levenberg-marquardt) methood, so it needs to be changed to differential evolution. These are my results: [![dielectric functions fitting by differential evolution method](http://juluribk.com/wp-content/uploads/2016/04/download.png){.aligncenter .size-full .wp-image-1614 width="839" height="337"}](http://juluribk.com/wp-content/uploads/2016/04/download.png) [You can download my ipython notebook at the bottom of this post.]{style="color: #ff0000;"} You can use this method for fitting other dielectric data, but have to tweak the bounds a little bit. If you use my implementation for your fitting your experimental data, please let me know or even better cite it as: Juluri B.K. "Differential Evolution algorithm to fit Metal Dielectric Functions" .
In the following code, I get the gold's experimental dielectric data from [Refractiveindex.info](http://refractiveindex.info/).
In $12$:
import urllib2 # To make a request and gather data, # Query the website response = urllib2.urlopen('http://refractiveindex.info/database/main/Au/Babar.yml') # Alternative Experimental data sources # response = urllib2.urlopen('http://refractiveindex.info/database/main/Au/McPeak.yml') # response = urllib2.urlopen('http://refractiveindex.info/database/main/Au/Lemarchand-3.96nm.yml') # Read the response. Response is in YAML format. Will use PyYAML parser, yaml_data = response.read() import yaml data = yaml.load(yaml_data) # Lets see what is inside the data, print"Keys in the data:" , data.keys() # Lets print the reference for this data, print "Reference : ", data['REFERENCES'] if 'COMMENTS' in data.keys(): print "Comments : ", data['COMMENTS'] # Data seems to be stored in 'DATA', if needed uncomment the following line #print(data['DATA'])
Keys in the data: ['REFERENCES', 'DATA'] Reference : S. Babar and J. H. Weaver. Optical constants of Cu, Ag, and Au revisited, Appl. Opt. 54, 477-481 (2015)
Data string contains real and imaginary part of the refractive index (not dielectric permitivitty). I clean up the real and imaginary part of the refractive data and convert them to dielectric permittivities
In $14$:
# Some clean up needed. The required data string holding wavelength, n and k, is in the form of dictionary item # with key 'data', which is inside a list\n", cleaned_nk = data['DATA'][0]['data'] # Use StringIO to convert String buffer as a file like class\n", import numpy as np from StringIO import StringIO # Read the data into as numpy vectors\n", wave_micron, n, k = np.genfromtxt(StringIO(cleaned_nk), unpack=True) wave_exp = wave_micron * 1000 # Convert wavelength from microns to nm eps_exp = (n + 1j* k)**2 # Convert n,k into dielectric function # Lets convert wavelengh in nm to eV and assign it to 'w' h = 4.135667516E-15 # plancks's constant in eV-sec c = 299792458E9 # speed of light in vacuum in nm/sec w_exp = h*c/wave_exp; # Convert wavleength in nm to eV
I define a function that will plot experimental dielectric data. It will also plot fitted models if the data is provided.
In $6$:
def plot_fit(w_exp,eps_exp, w_for_fit = None,eps_fit_result = None): ''' This functions plots the experimental dielectric function as function of wavelength. If fitted model data is available it will also plot it ''' %matplotlib inline import matplotlib.pyplot as plt from matplotlib.ticker import ScalarFormatter from matplotlib import rcParams rcParams.update({'font.size': 18}) # Increase the font size rcParams['mathtext.default']='regular' # make the mathtext the same size of normal text for better readbility # Plot the real part of the dielectric function fig,ax = plt.subplots(1,2,figsize = (12,5)) ax[0].scatter(w_exp, abs(eps_exp.real), marker = 'o',facecolor = 'none', edgecolor = 'g') ax[0].set_xlabel('Energy (eV)') ax[0].set_ylabel(r'|$\epsilon^\prime$|') # Plot the imaginary part of the dielectric function ax[1].scatter(w_exp, eps_exp.imag, marker = 'o',facecolor = 'none', edgecolor = 'g') ax[1].set_ylabel(r'$\epsilon ^ {\prime \prime}$') ax[1].set_xlabel('Energy (eV)') # If the fit data is available we will plot it with the actual plots if (w_for_fit is not None and eps_fit_result is not None): ax[0].plot(w_for_fit, abs(eps_fit_result.real),'-r') ax[1].plot(w_for_fit,eps_fit_result.imag,'-r') # Set the grid, log and axis formatter for axis in ax: axis.grid('on') axis.set_xscale('log') axis.set_yscale('log') axis.set_xlim([min(w_exp),max(w_exp)]) for x_yaxis in [axis.xaxis, axis.yaxis]: x_yaxis.set_major_formatter(ScalarFormatter()) fig.tight_layout()
Lets plot the experimental data to see the data. Data seems to be from 0.1ev to 10eV. There seems to be three transitions (are all of them interband?) above 1eV
In $7$:
plot_fit(w_exp,eps_exp)
In $15$:
In $9$:
In $17$: