Dear Meep users, I would like to calculate the scattering cross-section of a gold nanoparticle using Meep; with the method of hitting the nanoparticle with a plane wave and calculating the flux exiting a box that completely surrounds it. The calculation works reasonably if I model the nanoparticle with a Drude model: I haven't yet tried to refine the spatial grid, shorten the time step and so on but I get the scattering peak at the right position, the width also nicely compares with the width for pure dipolar scattering. The trouble starts when I want to do a Drude-Lorentz model with two additional Lorentz terms, in order to fit well the dielectric function of gold from 400 to 900 nm. The simulation then refuses to produce the expected scattering peak and the big rising shoulder at the blue side of the peak that is expected in scattering from a gold nanosphere (I see those in Mie calculations and also in a simple dipolar model).
I have tried to refine the time stepping, using a Courant factor of 0.1, but it has not helped at all. The curious thing is that on the opposite a two-dimensional model with the same parameters for the Drude-Lorentz model seems to work also nicely. The ctl file that I include here runs in about two hours on a rather fast machine (using four cores with mpi) - so obviously I would not ask anyone to run it ... but the point is, why does the model suddenly stop giving reasonable results? I had a suspicion that I may *still* have too long of a time-step using even the Courant factor of 0.1, but before trying this even longer simulation I place my question here. Besides, I have also used the nice ctl file at http://www.mail-archive.com/[email protected]/msg03527.html to verify my dielectric function, and that too seems correct. Where could the problem be? Thanks in advance for any answer, Giovanni ; This control file attempts to calculate the scattering from a dielectric nanosphere; it calculates the outgoing flux from a closed surface surrounding the nanosphere (define-param s_cell_x 4) ; size of cell in X direction (define-param s_cell_y 4) ; size of cell in Y direction (define-param s_cell_z 4) ; size of cell in Z direction (set! geometry-lattice (make lattice (size s_cell_x s_cell_y s_cell_z))) (set! Courant 0.5) ; define parameter for the presence of the sphere (define-param sphere_present? true) ; if true, the sphere is present, if false it is absent ; define parameters for the fluxes (define cx-1 -.75) (define cx0 0) (define cx1 .75) (define cy-1 -.75) (define cy0 0) (define cy1 .75) (define cz-1 -.75) (define cz0 0) (define cz1 .75) (define sx1 (* 2 cx1)) (define sx0 0) (define sy1 (* 2 cy1)) (define sy0 0) (define sz1 (* 2 cz1)) (define sz0 0) ; define gold (define gold (make dielectric (epsilon 4.98) (E-polarizations (make polarizability (omega 1e-10) (gamma 6.13e-4) (sigma (* .142 .142 1e20))) ; note that this is how one does arithmetical operations in the language Scheme. In the convention I use most often, * .14 .14 1e20 corresponds to .14 * .14 * 1e20 (make polarizability (omega 0.0592) (gamma 0.0211) (sigma 1.76)) (make polarizability (omega 0.0458) (gamma 0.0119) (sigma 0.952)) ) )) (define Drude_material (make dielectric (epsilon 1) (E-polarizations (make polarizability (omega 1e-10) (gamma 1e-3) (sigma (* .0377 .0377 3e20))) ; note that this is how one does arithmetical operations in the language Scheme. In the convention I use most often, * .14 .14 1e20 corresponds to .14 * .14 * 1e20 ) )) (set! geometry (if sphere_present? (list (make sphere (center 0 0 0) (radius .5) (material gold) ;(make dielectric (epsilon 1) ; (E-polarizations ; (make polarizability ; (omega 1e-10) (gamma 1e-3) (sigma (* .0377 .0377 3e20) ;))))) ; The parameters for the metal are chosen so that the resonance appears at a normalized frequency of about .00377 )) (list (make sphere (center 0 0 0) (radius .5) (material (make dielectric (epsilon 1))))) ; case with no sphere (the dielectric constant is set to 1 so that there is actually no sphere) ); end of if construct ); close parenthesis for !set geometry ; Set parameters for a gaussian source (define-param fcen .0377) ; pulse center frequency (define-param df .0377) ; pulse width (in frequency) (set! sources (list (make source (src (make gaussian-src (frequency fcen) (fwidth df) (cutoff 5))) (component Ez) (center -.9 0 0) (size 0 s_cell_y s_cell_z) ))) (set! symmetries (list (make mirror-sym (direction Y)))) (set! pml-layers (list (make pml (thickness 1.0)))) ;(set! pml-layers (list (make pml (thickness 3.0)))) (set! resolution 20) ;specify a plane through which Meep calculates a flux (define-param nfreq 100) ; number of frequencies at which to compute the flux (define scatt_left ; scattered flux through the left wall (add-flux fcen df nfreq (make flux-region (center cx-1 cy0 cz0) (size sx0 sy1 sz1) (weight -1.0)) ) ); close parenthesis for define scatt (define scatt_right ; scattered flux through the right wall (add-flux fcen df nfreq (make flux-region (center cx1 cy0 cz0) (size sx0 sy1 sz1) (weight 1.0)) ) ); close parenthesis for define scatt (define scatt_top ; scattered flux through the top wall (add-flux fcen df nfreq (make flux-region (center cx0 cy0 cz1) (size sx1 sy1 sz0) (weight 1.0)) ) ); close parenthesis for define scatt (define scatt_bottom ; scattered flux through the bottom wall (add-flux fcen df nfreq (make flux-region (center cx0 cy0 cz-1) (size sx1 sy1 sz0) (weight -1.0)) ) ); close parenthesis for define scatt (define scatt_front ; scattered flux through the bottom wall (add-flux fcen df nfreq (make flux-region (center cx0 cy-1 cz0) (size sx1 sy0 sz1) (weight -1.0)) ) ); close parenthesis for define scatt (define scatt_rear ; scattered flux through the bottom wall (add-flux fcen df nfreq (make flux-region (center cx0 cy1 cz0) (size sx1 sy0 sz1) (weight 1.0)) ) ); close parenthesis for define scatt (define inc ; incident flux (add-flux fcen df nfreq (make flux-region (center cx-1 cy0 cz0) (size 0 s_cell_y s_cell_z )) ) ); close parenthesis for define inc (use-output-directory) ;load the negated flux in case of no sphere so that we consider only the scattered flux (if sphere_present? (begin (load-minus-flux "direct-flux-left" scatt_left) (load-minus-flux "direct-flux-right" scatt_right) (load-minus-flux "direct-flux-top" scatt_top) (load-minus-flux "direct-flux-bottom" scatt_bottom) (load-minus-flux "direct-flux-front" scatt_front) (load-minus-flux "direct-flux-rear" scatt_rear)) ) (run-sources+ (stop-when-fields-decayed 50 Ez (vector3 .8 0 0) 1e-3) (at-beginning output-epsilon) ; (to-appended "ey" (at-every 0.6 output-efield-y)) ); close parenthesis for run-sources (if (not sphere_present?) (begin (save-flux "direct-flux-left" scatt_left) (save-flux "direct-flux-right" scatt_right) (save-flux "direct-flux-top" scatt_top) (save-flux "direct-flux-bottom" scatt_bottom) (save-flux "direct-flux-front" scatt_front) (save-flux "direct-flux-rear" scatt_rear)) ) (display-fluxes inc scatt_left scatt_top scatt_right scatt_bottom scatt_front scatt_rear) ---------------------------------------- Message sent by Cup Webmail (Roundcube) _______________________________________________ meep-discuss mailing list [email protected] http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss
