Hello Prof. and MPB users, I recently resumed my interest in MPB and began running simple calculations. I now look at 3D simple cubic grid of dielectric spheres. I run a ctl file where I calculate the bands, the fields and the power density. The title suggests that I might be overlooking something simple, but I will elaborate on my problem:
;define the cell geometry of spheres in cubic grid (define-param eps 11.9000000000000000000000000000000000000000000000000) (define-param r 0.2) (define-param dielectricon (make dielectric (epsilon eps) ) ) (set! default-material air) ;define the lattice vectors of cubic lattice (set! geometry-lattice (make lattice (basis-size 1 1 1) (basis1 1 0 0) (basis2 0 1 0) (basis3 0 0 1) ) ) (set! geometry (list (make sphere (center 0 0 0) (radius r) (material dielectricon) ) ) ) ;define the k-points (set! k-points (list (cartesian->reciprocal (vector3 0.2 0.3 -0.1)) ) ) ;define the simulation settings (set-param! resolution 32) (set-param! num-bands 1) (define-param top-bands 1) (set-param! tolerance 0.0000000001) ;disable default filename prefix in output files (set! filename-prefix false) ;special calculation functions to pass to run (define (output-nonbloch-efield which-band) (get-efield which-band) (cvector-field-nonbloch! cur-field) (output-field-to-file -1 "e.") ) (define (output-nonbloch-dfield which-band) (get-dfield which-band) (cvector-field-nonbloch! cur-field) (output-field-to-file -1 "d.") ) (define (output-nonbloch-hfield which-band) (get-hfield which-band) (cvector-field-nonbloch! cur-field) (output-field-to-file -1 "h.") ) (define (topbands-functions) (let topbands-loop ((n (- top-bands 1))) (output-nonbloch-efield (- num-bands n)) (output-nonbloch-dfield (- num-bands n)) (output-nonbloch-hfield (- num-bands n)) (output-tot-pwr (- num-bands n)) (output-hpwr (- num-bands n)) (output-dpwr (- num-bands n)) (if (> n 0) (topbands-loop (- n 1)) 0) ) ) (run topbands-functions) I than use matlab to read h5 files of the fields, compute the power density and write a new h5 file, then I use vis5d to visualize the the power densities. U = mpb.epsilon .* reshape(sum(abs(fields.E).^2), size(mpb.epsilon)) ... % fields.E is 4D matrix: 3xNxNxN, sum() will sum over the 3 components x,y,z squares Fx^2+Fy^2+Fz^2 + reshape(sum(abs(fields.H).^2), size(mpb.epsilon)); % epsilon is NxNxN matrix filename = [directory '\tot.rpwr.k' sprintf('%02d',k) '.b01.h5']; % or hpwr or dpwr in order to look at either the elect. or magn. energy Umpb = hdf5read(filename,'data'); file_new = [filename(1:end-3),'.U',filename(end-2:end)]; copyfile(filename, file_new); hdf5write(file_new,'/data',U); What I do not understand is as follows: 1. When I compare in vis5d (1) mpb output of the energy distribution and (2) my calculated energy distribution, I see that dpwr matches sum(abs(fields.E).^2) without multiplication in epsilon. It seems that output-nonbloch-efield gives D instead of E 2. when I check normalization of the fields using sum(abs(fields.F(:)).^2)*mpb.dV (dV being the spatial resolution of volume 1/N^3), I get that H is normalized to 1, E is normalized to 0.965 and D is normalized to 1.3333. According to the User Reference, D should also be normalized to 1, right? I checked my code again and again, so either there is a bug or I missed something too simple. Any help will be appreciated. Thanks in advance for the help.
_______________________________________________ mpb-discuss mailing list mpb-discuss@ab-initio.mit.edu http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/mpb-discuss