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
