In the past I've used the code below to do computations for wavevector
diagrams.
For example, to get a 2d wavevector diagram for the TM polarization in
a square or rectangular lattice, you would do:
(wavevector-diagram (kgrid 0 0.5 0 0.5 20 20) TM 0)
instead of (run-tm), then grep the output for lines beginning with
"kgrid:, 1," for the first band, "kgrid:, 2," for the second band, and
so on. You can then load this into Matlab or Octave or whatever and
do a contour plot etc.
That will only plot for 1/4 of the Brillouin zone, since I only
specified a k-point grid from 0 to 0.5 in the x and y directions. It
is usually easier to see what is going on if you plot for at least the
whole Brillouin zone, and sometimes for multiple copies of the
Brillouin zone. You can do that by specifying more k-points, of
course, but that is wasteful in the common case where you have some
rotational symmetry (so that the rest of the Brillouin zone is just
reflections of the upper-right quadrant). And, if you want multiple
copies of the Brillouin zone then of course you can just copy it. So
I wrote a little matlab function, replicate.m, to make copies of the
first quadrant of the Brillouin zone for me....see the Matlab code
below.
So, for example, I would read the kgrid data above into a variable f,
and then do:
contour(replicate(f, 3, 2))
to make 6 (3 x 2) copies of the Brillouin zone and then contour-plot it.
I would try something simple to start with, such as the isofrequency
diagram for a square lattice of dielectric rods in air, and compare to
the corresponding figure from the end of chapter 10 of our book
(readable online at ab-initio.mit.edu/book).
Steven
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Code to output the data needed for a wavevector diagram at a grid of
; k points, suitable for plotting with a contour-plot program.
(define (kgrid kx-min kx-max ky-min ky-max nkx nky)
(map (lambda (kx)
(interpolate nky (list (vector3 kx ky-min) (vector3 kx ky-
max))))
(interpolate nkx (list kx-min kx-max))))
; output a bunch of lines of the form:
; kgrid:, band#, kx, frequencies at kys...
; Frequencies above the light cone omega > c |k| / n-lightcone are
; excluded (multiplied by -1); set n-lightcone = 0 to disable this.
(define (wavevector-diagram kgrid parity n-lightcone)
(map (lambda (kylist)
(set! k-points kylist)
(run-parity parity true)
(map
(lambda (band)
(print "kgrid:, " band ", " (vector3-x (car kylist)))
(map (lambda (freqs k)
(print ", "
(* (if (and (positive? n-lightcone)
(> (list-ref freqs (- band 1))
(* n-lightcone
(vector3-norm
(reciprocal->cartesian
k)))))
-1
1)
(list-ref freqs (- band 1)))))
all-freqs k-points)
(print "\n"))
(arith-sequence 1 1 num-bands)))
kgrid))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function r = replicate(x, n, m)
x0 = x(:, 2:size(x,2)-1);
xu = [fliplr(x0), x];
xu0 = xu(2:size(xu,1)-1,:);
xs = [flipud(xu0); xu];
xsn = xs;
if mod(m,2) == 0
xsn = xs(:, size(x0,2)+1:size(xs,2));
end
for i = 2:m
xsn = [xsn, xs];
end
if mod(m,2) == 0
xsn = [xsn, xs(:, 1:size(x,2))];
end
r = xsn;
if mod(n,2) == 0
r = xsn(size(xu0,1)+1:size(xsn,1), :);
end
for j = 2:n
r = [r; xsn];
end
if mod(n,2) == 0
r = [r; xsn(1:size(xu,1), :)];
end
_______________________________________________
mpb-discuss mailing list
[email protected]
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/mpb-discuss
