Dear Steven, everybody,
I want to obtain the transmission/reflection spectra of a single gold
spherical nanoparticle located at the center of the simulation box.
Everything works fine except for the polarizability, which is not taken
into account. I use the Drude-Lorentz model described in Vial et al.,
PHYSICAL REVIEW B 71, 085416 (2005͒) with units adapted to my case. The
functions sigma_1 and sigma_2 here below are defined accordingly, and
the add-polarizability terms are called in the program. However, they
are not taken into account in the calculation, i.e. if I remove the
add-polarizability terms in the program, everything works the same way.
Owing to the fact that the same program written with the libctl
interface works perfectly, i.e. it reproduces adequately the expected
plasmon resonance peak(s), I guess something is wrong with my C++
implementation. Does anybody know why? The way the sigma_1 and sigma_2
functions are defined is not correct or is it due to the call of the
add-polarizability that is not fine?
Thanks in advance for your responses.
Renaud
#include <iostream>
#include <fstream>
#include "meep.hpp"
using namespace meep;
double epsgold=5.9673;
double resolution, D, fcen, fwidth;
const double pmlt=1.0, pad=3.0, padxy=1.0;
const double rh=0.3536;
const double sx=2*(padxy+pmlt+rh),sy=2*(padxy+pmlt+rh),sz=2*(pad+pmlt
+rh);
const double omega_d=D*sqrt(2.0)*7.045e-03; //plasma frequency of gold
in Drude model
const double sig1=1.0e+40, om1=omega_d*1.0e-20,
ga1=D*sqrt(2.0)*5.307e-05;
const double sig2=1.09, om2=D*sqrt(2.0)*2.170e-03,
ga2=D*sqrt(2.0)*3.495e-04;
complex<double> one(const vec &p) {
return 1.0;
}
double empty(const vec &p) {
return 1.0;
}
int read_(){
ifstream fin("parameters");
//PARAMETERS OF THE SIMULATION
fin >> resolution;
fin >> D;
fin >> fcen;
fin >> fwidth;
fin.close();
return 0;
}
double sigma_1(const vec &p) {
const grid_volume v = vol3d(sx,sy,sz,resolution);
vec xx = p - v.center();
if (sqrt((xx-vec(0,0,0))&(xx-vec(0,0,0))) < rh+1e-06) return sig1;
return 0.0;
}
double sigma_2(const vec &p) {
const grid_volume v = vol3d(sx,sy,sz,resolution);
vec xx = p - v.center();
if (sqrt((xx-vec(0,0,0))&(xx-vec(0,0,0))) < rh+1e-06) return sig2;
return 0.0;
}
double sphere(const vec &p) {
read_();
const grid_volume v = vol3d(sx,sy,sz,resolution);
vec xx = p - v.center();
if (sqrt((xx-vec(0,0,0))&(xx-vec(0,0,0))) < rh+1e-06) return epsgold;
return 1.0;
}
void excitefilm(double a, component c) {
grid_volume v = vol3d(sx,sy,sz,a);
structure s0(v,empty,pml(pmlt));
structure s(v,sphere,pml(pmlt));
s.add_polarizability(sigma_1,om1,ga1);
s.add_polarizability(sigma_2,om2,ga2);
fields f0(&s0);
fields f(&s);
f0.use_real_fields();
f.use_real_fields();
f0.use_bloch(vec(0,0,0));
f.use_bloch(vec(0,0,0));
f.output_hdf5(Dielectric, v.surroundings());
int nfreq = 500;
gaussian_src_time src(fcen, fwidth);
volume src_plane(vec(pmlt+0.5,pmlt
+0.5,sz-pmlt-0.5),vec(sx-pmlt-0.5,sy-pmlt-0.5,sz-pmlt-0.5));
volume Rflux_plane(vec(pmlt+0.5,pmlt
+0.5,sz-pmlt-0.5),vec(sx-pmlt-0.5,sy-pmlt-0.5,sz-pmlt-0.5));
volume Tflux_plane(vec(pmlt+0.5,pmlt+0.5,pmlt
+0.5),vec(sx-pmlt-0.5,sy-pmlt-0.5,pmlt+0.5));
f0.add_volume_source(c,src,src_plane,one,1.0);
f.add_volume_source(c,src,src_plane,one,1.0);
master_printf("volume sources added...\n");
dft_flux ft0 = f0.add_dft_flux_plane(Tflux_plane,fcen-0.5*fwidth,fcen
+0.5*fwidth,nfreq);
dft_flux fr0 = f0.add_dft_flux_plane(Rflux_plane,fcen-0.5*fwidth,fcen
+0.5*fwidth,nfreq);
dft_flux ft = f.add_dft_flux_plane(Tflux_plane,fcen-0.5*fwidth,fcen
+0.5*fwidth,nfreq);
dft_flux fr = f.add_dft_flux_plane(Rflux_plane,fcen-0.5*fwidth,fcen
+0.5*fwidth,nfreq);
master_printf("simulating empty structure...\n");
while (f0.time() < f0.last_source_time() + 100) f0.step();
fr -= fr0;
master_printf("simulating sphere...\n");
while (f.time() < f.last_source_time() + 100) {
f.step();
if (int(f.time())%200 == 0) f.output_hdf5(Ex, v.surroundings());
}
double *flux = ft.flux();
double *flux0 = ft0.flux();
double *T;
double *R;
T = new double[nfreq];
R = new double[nfreq];
for (int i=0; i<nfreq; ++i)
T[i] = flux[i] / flux0[i];
delete[] flux;
flux = fr.flux();
for (int i=0; i<nfreq; ++i)
R[i] = -flux[i] / flux0[i];
delete [] flux;
delete [] flux0;
double dfreq = fwidth / (nfreq-1);
master_printf("tranmission:, freq, R[i], T[i]\n");
for (int i=0; i<nfreq; ++i)
master_printf("transmission:, %f, %f, %f\n",fcen-0.5*fwidth
+i*dfreq,R[i],T[i]);
delete [] R;
delete [] T;
}
int main(int argc, char **argv) {
initialize mpi(argc, argv);
master_printf("beginning reflection/tranmission spectra of gold
sphere ...\n");
excitefilm(80, Ex);
master_printf("finished.\n");
return 0;
}
--
Dr. Renaud Vallée
Centre de Recherche Paul Pascal
115, avenue du docteur Albert Schweitzer
33600 Pessac
France
Tél: +33556845612
_______________________________________________
meep-discuss mailing list
[email protected]
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss
