-
Notifications
You must be signed in to change notification settings - Fork 85
/
user-ref.html
1725 lines (1390 loc) · 74.9 KB
/
user-ref.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 3.2//EN">
<HTML><HEAD>
<TITLE>User Reference</TITLE>
<LINK rel="Contents" href="index.html">
<LINK rel="Copyright" href="license.html">
<LINK rel="Start" href="index.html">
</HEAD>
<BODY TEXT="#000000" BGCOLOR="#FFFFFF">
Go to the <a href="developer.html">next</a>, <a href="analysis-tutorial.html">previous</a>, or <a href="index.html">main</a> section.
<hr>
<h1>User Reference</h1>
<p>Here, we document the features exposed to the user by the MIT
Photonic-Bands package. We do not document the Scheme language or the
functions provided by libctl (see also the <a
href="http://ab-initio.mit.edu/libctl/doc/user-ref.html">user
reference</a> section of the <a
href="http://ab-initio.mit.edu/libctl/doc/">libctl manual</a>).
<h2><a name="input-vars">Input Variables</a></h2>
<p>These are global variables that you can set to control various
parameters of the Photonic-Bands computation. They are also listed
(along with their current values) by the <code>(help)</code> command.
In brackets after each variable is the type of value that it should
hold. (The classes, complex datatypes like
<code>geometric-object</code>, are described in a later subsection.
The basic datatypes, like <code>integer</code>, <code>boolean</code>,
<code>cnumber</code>, and <code>vector3</code>, are defined by libctl.)
<dl>
<dt><code>geometry</code> [list of <code>geometric-object</code> class]
<dd>Specifies the geometric objects making up the structure being
simulated. When objects overlap, later objects in the list take
precedence. Defaults to no objects (empty list).
<p><dt><code>default-material</code> [<code>material-type</code> class]
<dd>Holds the default material that is used for points not in any
object of the geometry list. Defaults to air (epsilon of 1). See
also <code>epsilon-input-file</code>, below.
<p><dt><code>ensure-periodicity</code> [<code>boolean</code>]
<dd>If true (the default), then geometric objects are treated as if
they were shifted by all possible lattice vectors; i.e. they are made
periodic in the lattice.
<p><dt><code>geometry-lattice</code> [<code>lattice</code> class]
<dd>Specifies the basis vectors and lattice size of the computational
cell (which is centered on the origin of the coordinate system).
These vectors form the basis for all other 3-vectors in the geometry,
and the lattice size determines the size of the primitive cell. If
any dimension of the lattice size is the special value
<code>no-size</code>, then the dimension of the lattice is reduced
(i.e. it becomes two- or one-dimensional). (That is, the dielectric
function becomes two-dimensional; it is still, in principle, a three
dimensional system, and the k-point vectors can be three-dimensional.)
Generally, you should make any <code>no-size</code> dimension(s)
perpendicular to the others. Defaults to the orthogonal x-y-z vectors
of unit length (i.e. a square/cubic lattice).
<p><dt><code>resolution</code> [<code>number</code> or <code>vector3</code>]
<dd>Specifies the computational grid resolution, in pixels per lattice
unit (a lattice unit is one basis vector in a given direction). If
<code>resolution</code> is a <code>vector3</code>, then specifies a
different resolution for each direction; otherwise the resolution is
uniform. (The grid size is then the product of the lattice size and
the resolution, rounded up to the next positive integer.) Defaults to
<code>10</code>.
<p><dt><code>grid-size</code> [<code>vector3</code>]
<dd>Specifies the size of the discrete computational grid along each
of the lattice directions. <em>Deprecated:</em> the preferred method
is to use the <code>resolution</code> variable, above, in which case
the <code>grid-size</code> defaults to <code>false</code>. To get the
grid size you should instead use the <code>(get-grid-size)</code>
function.
<p><dt><code>mesh-size</code> [<code>integer</code>]
<dd>At each grid point, the dielectric constant is averaged over a
"mesh" of points to find an effective dielectric tensor. This mesh is
a grid with <code>mesh-size</code> points on a side. Defaults to 3.
Increasing <code>mesh-size</code> makes the the average dielectric
constant sensitive to smaller structural variations without increasing
the grid size, but also means that computing the dielectric function
will take longer. (Using a <code>mesh-size</code> of <code>1</code>
turns <em>off</em> dielectric function averaging and the creation of
an effective dielectric tensor at interfaces. Be sure you know what
you're doing before you worsen convergence in this way.)
<p><dt><code>dimensions</code> [<code>integer</code>]
<dd>Explicitly specifies the dimensionality of the simulation; if the
value is less than 3, the sizes of the extra dimensions in
<code>grid-size</code> are ignored (assumed to be one). Defaults to 3.
<em>Deprecated:</em> the preferred method is to set
<code>geometry-lattice</code> to have size </code>no-size</code> in
any unwanted dimensions.
<p><dt><code><a name="k-points">k-points</a></code> [list of
<code>vector3</code>]
<dd>List of Bloch wavevectors to compute the bands at, expressed in
the basis of the reciprocal lattice vectors. The reciprocal lattice
vectors are defined as follows: Given the lattice vectors
<code>R<sub>i</sub></code> (<em>not</em> the basis vectors), the
reciprocal lattice vector <code>G<sub>j</sub></code> satisfies
<code>R<sub>i</sub> * G<sub>j</sub> = 2*Pi*delta<sub>i,j</sub></code>,
where <code>delta<sub>i,j</sub></code> is the Kronecker delta (1 for
<code>i = j</code> and 0 otherwise). (<code>R<sub>i</sub></code> for
any <code>no-size</code> dimensions is taken to be the corresponding
basis vector.) Normally, the wavevectors should be in the first
Brillouin zone (<a href="#first-brillouin-zone">see below</a>).
<code>k-points</code> defaults to none (empty list).
<p><dt><code>num-bands</code> [<code>integer</code>]
<dd>Number of bands (eigenvectors) to compute at each k
point. Defaults to 1.
<p><dt><code>target-freq</code> [<code>number</code>]
<dd>If zero, the lowest-frequency <code>num-bands</code> states are
solved for at each k point (ordinary eigenproblem). If non-zero,
solve for the <code>num-bands</code> states whose frequencies are have
the smallest absolute difference with <code>target-freq</code>
(special, "targeted" eigenproblem). Beware that the targeted solver
converges more slowly than the ordinary eigensolver and may require a
lower <code>tolerance</code> to get reliable results. Defaults to 0.
<p><dt><code>tolerance</code> [<code>number</code>]
<dd>Specifies when convergence of the eigensolver is judged to have
been reached (when the eigenvalues have a fractional change less than
<code>tolerance</code> between iterations). Defaults to 1.0e-7.
<p><dt><code>filename-prefix</code> [<code>string</code>]
<dd>A string prepended to all output filenames. Defaults to
<code>""</code> (no prefix).
<p><dt><code>epsilon-input-file</code> [<code>string</code>]
<dd>If this string is not <code>""</code> (the default), then it
should be the name of an HDF5 file whose first/only dataset defines a
dielectric function (over some discrete grid). This dielectric
function is then used in place of <code>default-material</code>
(<i>i.e.</i> where there are no <code>geometry</code> objects). The
grid of the epsilon file dataset need not match
<code>grid-size</code>; it is scaled and/or linearly interpolated as
needed. The lattice vectors for the epsilon file are assumed to be
the same as <code>geometry-lattice</code>. [ Note that, even if the
grid sizes match and there are no geometric objects, the dielectric
function used by MPB will not be exactly the dielectric function of
the epsilon file, unless you also set <code>mesh-size</code> to
<code>1</code> (see above). ]
<p><dt><code>eigensolver-block-size</code> [<code>integer</code>]
<dd>The eigensolver uses a "block" algorithm, which means that it
solves for several bands simultaneously at each k-point.
<code>eigensolver-block-size</code> specifies this number of bands to
solve for at a time; if it is zero or >= <code>num-bands</code>,
then all the bands are solved for at once. If
<code>eigensolver-block-size</code> is a negative number, -<i>n</i>,
then MPB will try to use nearly-equal block-sizes close to <i>n</i>.
Making the block size a small number can reduce the memory
requirements of MPB, but block sizes > 1 are usually more efficient
(there is typically some optimum size for any given problem).
Defaults to -11 (i.e. solve for around 11 bands at a time).
<p><dt><code>simple-preconditioner?</code> [<code>boolean</code>]
<dd>Whether or not to use a simplified preconditioner; defaults to
<code>false</code> (this is fastest most of the time). (Turning this
on increases the number of iterations, but decreases the time for each
iteration.)
<p><dt><code>deterministic?</code> [<code>boolean</code>]
<dd>Since the fields are initialized to random values at the start of
each run, there are normally slight differences in the number of
iterations, etcetera, between runs. Setting
<code>deterministic?</code> to <code>true</code> makes things
deterministic (the default is <code>false</code>).
<p><dt><code>eigensolver-flags</code> [<code>integer</code>]
<dd>This variable is undocumented and reserved for use by Jedi Masters only.
</dl>
<h2><a name="predef-vars">Predefined Variables</a></h2>
<p>Variables predefined for your convenience and amusement.
<dl>
<dt><code>air</code>, <code>vacuum</code> [<code>material-type</code> class]
<dd>Two aliases for a predefined material type with a dielectric
constant of 1.
<dt><code>nothing</code> [<code>material-type</code> class]
<dd>A material that, effectively, punches a hole through other
objects to the background (<code>default-material</code> or
<code>epsilon-input-file</code>).
<p><dt><code>infinity</code> [<code>number</code>]
<dd>A big number (1.0e20) to use for "infinite" dimensions of objects.
</dl>
<h2><a name="output-vars">Output Variables</a></h2>
<p>Global variables whose values are set upon completion of the
eigensolver.
<dl>
<dt><code>freqs</code> [list of <code>number</code>]
<dd>A list of the frequencies of each band computed for the last k
point. Guaranteed to be sorted in increasing order. The frequency of
band <code>b</code> can be retrieved via <code>(list-ref freqs (- b
1))</code>.
<p><dt><code>iterations</code> [<code>integer</code>]
<dd>The number of iterations required for convergence of the last k point.
<p><dt><code>parity</code> [<code>string</code>]
<dd>A string describing the current required parity/polarization
(<code>"te"</code>, <code>"zeven"</code>, etcetera, or <code>""</code>
for none), useful for prefixing output lines for grepping.
</dl>
<p>Yet more global variables are set by the <code>run</code> function
(and its variants), for use after <code>run</code> completes or by a
band function (which is called for each band during the execution of
<code>run</code>.
<dl>
<dt><code>current-k</code> [<code>vector3</code>]
<dd>The k point (from the <code>k-points</code> list) most recently
solved.
<p><dt><code>gap-list</code> [list of <code>(<i>percent freq-min freq-max</i>)</code> lists]
<dd>This is a list of the gaps found by the eigensolver, and is set by
the <code>run</code> functions when two or more k-points are solved.
(It is the empty list if no gaps are found.)
<p><dt><code>band-range-data</code> [list of <code>((<i>min . kpoint</i>) . (<i>max . kpoint</i>))</code> pairs (of pairs)]
<dd>For each band, this list contains the minimum and maximum
frequencies of the band, and the associated k points where the extrema
are achieved. Note that the bands are defined by sorting the
frequencies in increasing order, so this can be confused if two bands
cross.
</dl>
<h2><a name="classes">Classes</a></h2>
<p>Classes are complex datatypes with various "properties" which may
have default values. Classes can be "subclasses" of other classes;
subclasses inherit all the properties of their superclass, and can be
used any place the superclass is expected. An object of a class is
constructed with:
<pre>
(make <i>class</i> (<i>prop1 val1</i>) (<i>prop2 val2</i>) ...)
</pre>
<p>See also the <a href="http://ab-initio.mit.edu/libctl/doc/">libctl
manual</a>.
<p>MIT Photonic-Bands defines several types of classes, the most
numerous of which are the various geometric object classes. You can
also get a list of the available classes, along with their property
types and default values, at runtime with the <code>(help)</code>
command.
<h3><a name="lattice">lattice</a></h3>
<p>The lattice class is normally used only for the
<code>geometry-lattice</code> variable, and specifies the three
lattice directions of the crystal and the lengths of the corresponding
lattice vectors.
<dl>
<dt><code>lattice</code>
<dd>Properties:
<dl>
<dt><code>basis1</code>, <code>basis2</code>, <code>basis3</code> [<code>vector3</code>]
<dd>The three lattice directions of the crystal, specified in the
cartesian basis. The lengths of these vectors are ignored--only their
directions matter. The lengths are determined by the
<code>basis-size</code> property, below. These vectors are then used
as a basis for all other 3-vectors in the ctl file. They default to
the x, y, and z directions, respectively.
<dt><code>basis-size</code> [<code>vector3</code>]
<dd>The components of <code>basis-size</code> are the lengths of the
three basis vectors, respectively. They default ot unit lengths.
<dt><code>size</code> [<code>vector3</code>]
<dd>The size of the lattice (i.e. the length of the lattice vectors
<code>R<sub>i</sub></code>, in which the crystal is periodic) in units
of the basis vectors. Thus, the actual lengths of the lattice vectors
are given by the components of <code>size</code> multiplied by the
components of <code>basis-size</code>. (Alternatively, you can think
of <code>size</code> as the vector between opposite corners of the
primitive cell, specified in the lattice basis.) Defaults to unit
lengths.
<p>If any dimension has the special size <code>no-size</code>, then
the dimensionality of the problem is reduced by one; strictly
speaking, the dielectric function is taken to be uniform along that
dimension. (In this case, the <code>no-size</code> dimension should
generally be orthogonal to the other dimensions.)
</dl>
</dl>
<h3><a name="material-type">material-type</a></h3>
<p>This class is used to specify the materials that geometric objects
are made of. Currently, there are three subclasses,
<code>dielectric</code>, <code>dielectric-anisotropic</code>, and
<code>material-function</code>.
<dl>
<dt><code>dielectric</code>
<dd>A uniform, isotropic, linear dielectric material, with one property:
<dl>
<dt><code>epsilon</code> [<code>number</code>]
<dd>The dielectric constant (must be positive). No default value. You
can also use <code>(index <i>n</i>)</code> as a synonym for
<code>(epsilon (* <i>n n</i>))</code>.
</dl>
<p><dt><code><a name="dielectric-anisotropic">dielectric-anisotropic</a></code>
<dd>A uniform, possibly anisotropic, linear dielectric material. For
this material type, you specify the (real-symmetric, or possibly
complex-hermitian) dielectric tensor (relative to the cartesian xyz
axes):
<pre>
a u v
epsilon = ( u* b w )
v* w* c
</pre>
<p>This allows your dielectric to have different dielectric constants
for fields polarized in different directions. The epsilon tensor must
be positive-definite (have all positive eigenvalues); if it is not,
MPB exits with an error. (This does <em>not</em> imply that all of
the entries of the epsilon matrix need be positive.)
<p>The components of the tensor are specified via three properties:
<dl>
<dt><code>epsilon-diag</code> [<code>vector3</code>]
<dd>The diagonal elements (a b c) of the dielectric tensor. No default value.
<dt><code>epsilon-offdiag</code> [<code>cvector3</code>]
<dd>The off-diagonal elements (u v w) of the dielectric tensor.
Defaults to zero. This is a <code>cvector3</code>, which simply
means that the components may be complex numbers
(e.g. <code>3+0.1i</code>). If non-zero imaginary parts are
specified, then the dielectric tensor is complex-hermitian. This is
only supported when MPB is configured with the
<code>--with-hermitian-eps</code> flag. This is not dissipative (the
eigenvalues of epsilon are real), but rather breaks time-reversal
symmetry, corresponding to a gyrotropic (magneto-optic) material.
Note that <a href="#inv-symmetry">inversion symmetry</a> may not mean
what you expect for complex-hermitian epsilon, so be cautious about
using <code>mpbi</code> in this case.
<dt><code>epsilon-offdiag-imag</code> [<code>vector3</code>]
<dd><em>Deprecated:</em> The imaginary parts of the off-diagonal
elements (u v w) of the dielectric tensor; defaults to zero. Setting
the imaginary parts directly by specifying complex numbers in
<code>epsilon-offdiag</code> is preferred.
</dl>
<p>For example, a material with a dielectric constant of 3.0 for TE
fields (polarized in the xy plane) and 5.0 for TM fields (polarized in
the z direction) would be specified via <code>(make
(dielectric-anisotropic (epsilon-diag 3 3 5)))</code>. Please <a
href="#tetm-aniso">be aware</a> that not all 2d anisotropic dielectric
structures will have TE and TM modes, however.
<p><dt><code>material-function</code>
<dd>This material type allows you to specify the material as an
arbitrary function of position. (For an example of this, see the
<code>bragg-sine.ctl</code> file in the <code>examples/</code>
directory.) It has one property:
<dl>
<dt><code>material-func</code> [<code>function</code>]
<dd>A function of one argument, the position <code>vector3</code> (in
lattice coordinates), that returns the material at that point.
</dl>
<p>Note that the function you supply can return <em>any</em> material;
wild and crazy users could even return another
<code>material-function</code> object (which would then have its
function invoked in turn).
</dl>
<p>Normally, the dielectric constant is required to be positive (or
positive-definite, for a tensor). However, MPB does have a somewhat
experimental feature allowing negative dielectrics (e.g. in a plasma).
To use it, call the function <code>(allow-negative-epsilon)</code>
before <code>(run)</code>. In this case, it will output the (real)
frequency <em>squared</em> in place of the (possibly imaginary)
frequencies. (Convergence will be somewhat slower because the
eigenoperator is not positive definite.)
<h3><a name="geometric-object">geometric-object</a></h3>
<p>This class, and its descendants, are used to specify the solid
geometric objects that form the dielectric structure being simulated.
The base class is:
<dl>
<dt><code>geometric-object</code>
<dd>Properties:
<dl>
<dt><code>material</code> [<code>material-type</code> class]
<dd>The material that the object is made of (usually some sort of
dielectric). No default value (must be specified).
<dt><code>center</code> [<code>vector3</code>]
<dd>Center point of the object. No default value.
</dl>
</dl>
<p>One normally does not create objects of type
<code>geometric-object</code> directly, however; instead, you use one
of the following subclasses. Recall that subclasses inherit the
properties of their superclass, so these subclasses automatically have
the <code>material</code> and <code>center</code> properties (which
must be specified, since they have no default values).
<p>Recall that all 3-vectors, including the center of an object, its
axes, and so on, are specified in the basis of the normalized lattice
vectors normalized to <code>basis-size</code>. Note also that
3-vector properties can be specified by either <code>(<i>property</i>
(vector3 <i>x y z</i>))</code> or, equivalently,
<code>(<i>property</i> <i>x y z</i>)</code>.
<p>In a two-dimensional calculation, only the intersections of the
objects with the x-y plane are considered.
<dl>
<dt><code>sphere</code>
<dd>A sphere. Properties:
<dl>
<dt><code>radius</code> [<code>number</code>]
<dd>Radius of the sphere. No default value.
</dl>
<p><dt><code>cylinder</code>
<dd>A cylinder, with circular cross-section and finite height. Properties:
<dl>
<dt><code>radius</code> [<code>number</code>]
<dd>Radius of the cylinder's cross-section. No default value.
<dt><code>height</code> [<code>number</code>]
<dd>Length of the cylinder along its axis. No default value.
<dt><code>axis</code> [<code>vector3</code>]
<dd>Direction of the cylinder's axis; the length of this vector is
ignored. Defaults to point parallel to the z axis.
</dl>
<p><dt><code>cone</code>
<dd>A cone, or possibly a truncated cone. This is actually a subclass
of <code>cylinder</code>, and inherits all of the same properties,
with one additional property. The radius of the base of the cone is
given by the <code>radius</code> property inherited from
<code>cylinder</code>, while the radius of the tip is given by the new
property:
<dl>
<dt><code>radius2</code> [<code>number</code>]
<dd>Radius of the tip of the cone (i.e. the end of the cone pointed to
by the <code>axis</code> vector). Defaults to zero (a "sharp" cone).
</dl>
<p><dt><code>block</code>
<dd>A parallelepiped (i.e., a brick, possibly with non-orthogonal
axes). Properties:
<dl>
<dt><code>size</code> [<code>vector3</code>]
<dd>The lengths of the block edges along each of its three axes. Not
really a 3-vector (at least, not in the lattice basis), but it has
three components, each of which should be nonzero. No default value.
<dt><code>e1</code>, <code>e2</code>, <code>e3</code> [<code>vector3</code>]
<dd>The directions of the axes of the block; the lengths of these
vectors are ignored. Must be linearly independent. They default to
the three lattice directions.
</dl>
<p><dt><code>ellipsoid</code>
<dd>An ellipsoid. This is actually a subclass of <code>block</code>,
and inherits all the same properties, but defines an ellipsoid
inscribed inside the block.
</dl>
<p>Here are some examples of geometric objects created using the above
classes, assuming that the lattice directions (the basis) are just the
ordinary unit axes, and <code>m</code> is some material we have
defined:
<pre>
; A cylinder of infinite radius and height 0.25 pointing along the x axis,
; centered at the origin:
(make cylinder (center 0 0 0) (material m)
(radius infinity) (height 0.25) (axis 1 0 0))
; An ellipsoid with its long axis pointing along (1,1,1), centered on
; the origin (the other two axes are orthogonal and have equal
; semi-axis lengths):
(make ellipsoid (center 0 0 0) (material m)
(size 0.8 0.2 0.2)
(e1 1 1 1)
(e2 0 1 -1)
(e3 -2 1 1))
; A unit cube of material m with a spherical air hole of radius 0.2 at
; its center, the whole thing centered at (1,2,3):
(set! geometry (list
(make block (center 1 2 3) (material m) (size 1 1 1))
(make sphere (center 1 2 3) (material air) (radius 0.2))))
</pre>
<h2><a name="functions">Functions</a></h2>
<p>Here, we describe the functions that are defined by the
Photonic-Bands package. There are many types of functions defined,
ranging from utility functions for duplicating geometric objects to
run functions that start the computation.
<p>See also the <a
href="http://ab-initio.mit.edu/libctl/doc/user-ref.html">reference
section</a> of the libctl manual, which describes a number of useful
functions defined by libctl.
<h3><a name="geom-utils">Geometry utilities</a></h3>
<p>Some utility functions are provided to help you manipulate
geometric objects:
<dl>
<dt><code>(shift-geometric-object <i>obj shift-vector</i>)</code>
<dd>Translate <code>obj</code> by the 3-vector <code>shift-vector</code>.
<p><dt><code>(geometric-object-duplicates <i>shift-vector min-multiple max-multiple obj</i>)</code>
<dd>Return a list of duplicates of <code>obj</code>, shifted by
various multiples of <code>shift-vector</code> (from
<code>min-multiple</code> to <code>max-multiple</code>, inclusive, in
steps of 1).
<p><dt><code>(geometric-objects-duplicates <i>shift-vector min-multiple max-multiple obj-list</i>)</code>
<dd>Same as <code>geometric-object-duplicates</code>, except operates
on a list of objects, <code>obj-list</code>. If <i>A</i> appears
before <i>B</i> in the input list, then all the duplicates of <i>A</i>
appear before all the duplicates of <i>B</i> in the output list.
<p><dt><code>(geometric-objects-lattice-duplicates <i>obj-list [ ux uy uz ]</i>)</code>
<dd>Duplicates the objects in <code>obj-list</code> by multiples of
the lattice basis vectors, making all possible shifts of the
"primitive cell" (see below) that fit inside the lattice cell. (This
is useful for supercell calculations; see the <a
href="user-tutorial.html">tutorial</a>.) The primitive cell to
duplicate is <code>ux</code> by <code>uy</code> by <code>uz</code>, in
units of the basis vectors. These three parameters are optional; any
that you do not specify are assumed to be <code>1</code>.
<p><dt><code>(point-in-object? <i>point obj</i>)</code>
<dd>Returns whether or not the given 3-vector <code>point</code> is
inside the geometric object <code>obj</code>.
<p><dt><code>(point-in-periodic-object? <i>point obj</i>)</code>
<dd>As <code>point-in-object?</code>, but also checks translations of
the given object by the lattice vectors.
<p><dt><code>(display-geometric-object-info <i>indent-by obj</i>)</code>
<dd>Outputs some information about the given <code>obj</code>,
indented by <code>indent-by</code> spaces.
</dl>
<h3><a name="coordconvert">Coordinate conversion functions</a></h3>
<p>The following functions allow you to easily convert back and forth
between the lattice, cartesian, and reciprocal bases. (See also the
<a href="user-tutorial.html#units">note on units</a> in the tutorial.)
<dl>
<dt><code>(lattice->cartesian <i>x</i>)</code>, <code>(cartesian->lattice <i>x</i>)</code>
<dd>Convert <code><i>x</i></code> between the lattice basis (the basis
of the lattice vectors normalized to <code>basis-size</code>) and the
ordinary cartesian basis, where <code><i>x</i></code> is either a
<code>vector3</code> or a <code>matrix3x3</code>, returning the
transformed vector/matrix. In the case of a matrix argument, the
matrix is treated as an operator on vectors in the given basis, and is
transformed into the same operator on vectors in the new basis.
<p><dt><code>(reciprocal->cartesian <i>x</i>)</code>, <code>(cartesian->reciprocal <i>x</i>)</code>
<dd>Like the above, except that they convert to/from reciprocal space
(the basis of the reciprocal lattice vectors). Also, the cartesian
vectors output/input are in units of 2 Pi.
<p><dt><code>(reciprocal->lattice <i>x</i>)</code>, <code>(lattice->reciprocal <i>x</i>)</code>
<dd>Convert between the reciprocal and lattice bases, where the
conversion again leaves out the factor of 2 Pi (i.e. the lattice-basis
vectors are assumed to be in units of 2 Pi).
</dl>
<p>Also, a couple of rotation functions are defined, for convenience,
so that you don't have to explicitly convert to cartesian coordinates
in order to use libctl's <code>rotate-vector3</code> function (see the
<a href="http://ab-initio.mit.edu/libctl/doc/user-ref.html">libctl
reference</a>):
<dl>
<dt><code>(rotate-lattice-vector3 <i>axis theta v</i>)</code>, <code>(rotate-reciprocal-vector3 <i>axis theta v</i>)</code>
<dd>Like <code>rotate-vector3</code> , except that
<code><i>axis</i></code> and <code><i>v</i></code> are specified in
the lattice/reciprocal bases.
</dl>
<p><a name="first-brillouin-zone"></a>Usually, k-points are specified
in the first Brillouin zone, but sometimes it is convenient to specify
an arbitrary k-point. However, the accuracy of MPB degrades as you
move farther from the first Brillouin zone (due to the choice of a
fixed planewave set for a basis). This is easily fixed: simply
transform the k-point to a corresponding point in the first Brillouin
zone, and a completely equivalent solution (identical frequency,
fields, etcetera) is obtained with maximum accuracy. The following
function accomplishes this:
<dl>
<dt><code>(first-brillouin-zone <i>k</i>)</code>
<dd>Given a k-point <code><i>k</i></code> (in the basis of the
reciprocal lattice vectors, as usual), return an equivalent point in
the first Brillouin zone of the current lattice
(<code>geometry-lattice</code>).
</dl>
<p>Note that <code>first-brillouin-zone</code> can be applied to the
entire <code>k-points</code> list with the Scheme expression:
<code>(map first-brillouin-zone k-points)</code>.
<h3><a name="run">Run functions</a></h3>
<p>These are functions to help you run and control the simulation.
The ones you will most commonly use are the <code>run</code> function
and its variants. The syntax of these functions, and one lower-level
function, is:
<dl>
<dt><code>(run <i>band-func ...</i>)</code> <dd>This runs the
simulation described by the input parameters (see above), with no
constraints on the polarization of the solution. That is, it reads
the input parameters, initializes the simulation, and solves for the
requested eigenstates of each k-point. The dielectric function is
outputted to "<code>epsilon.h5</code>" before any eigenstates are
computed. <code>run</code> takes as arguments zero or more "band
functions" <code>band-func</code>. A band function should be a
function of one integer argument, the band index, so that
<code>(band-func which-band)</code> performs some operation on the
band <code>which-band</code> (e.g. outputting fields). After every
k-point, each band function is called for the indices of all the bands
that were computed. Alternatively, a band function may be a "thunk"
(function of zero arguments), in which case <code>(band-func)</code>
is called exactly once per k-point.
<p><dt><code>(run-zeven <i>band-func ...</i>), (run-zodd <i>band-func ...</i>)</code>
<dd>These are the same as the <code>run</code> function except that
they constrain their solutions to have even and odd symmetry with
respect to the z=0 plane. You should use these functions
<em>only</em> for structures that are symmetric through the z=0 mirror
plane, where the third basis vector is in the z direction (0,0,1) and
is orthogonal to the other two basis vectors, and when the k vectors
are in the xy plane. Under these conditions, the eigenmodes always
have either even or odd symmetry. In two dimensions, even/odd
parities are equivalent to TE/TM polarizations, respectively (and are
often strongly analogous even in 3d). Such a symmetry classification
is useful for structures such as waveguides and photonic-crystal slabs.
(For example, see the paper by S. G. Johnson <i>et al.</i>, "Guided
modes in photonic crystal slabs," <i>PRB</i> <b>60</b>, 5751, August
1999.)
<p><dt><code>(run-te <i>band-func ...</i>), (run-tm <i>band-func ...</i>)</code>
<dd>These are the same as the <code>run</code> function except that
they constrain their solutions to be TE- and TM-polarized,
respectively, in two dimensions. The TE and TM polarizations are
defined has having electric and magnetic fields in the xy plane,
respectively. Equivalently, the H/E field of TE/TM light has only a z
component (making it easier to visualize).
<p>These functions are actually equivalent to calling
<code>run-zeven</code> and <code>run-zodd</code>, respectively.
<p><a name="tetm-aniso">Note</a> that for the modes to be segregated
into TE and TM polarizations, the dielectric function must have mirror
symmetry for reflections through the xy plane. If you use <a
href="#dielectric-anisotropic">anisotropic dielectrics</a>, you should
be aware that they break this symmetry if the z direction is not one
of the principle axes. If you use <code>run-te</code> or
<code>run-tm</code> in such a case of broken symmetry, MPB will exit
with an error.
<p><dt><code>(run-yeven <i>band-func ...</i>), (run-yodd <i>band-func ...</i>)</code>
<dd>These functions are analogous to <code>run-zeven</code> and
<code>run-zodd</code>, except that they constrain their solutions to
have even and odd symmetry with respect to the y=0 plane. You should
use these functions <em>only</em> for structures that are symmetric
through the y=0 mirror plane, where the second basis vector is in the
y direction (0,1,0) and is orthogonal to the other two basis vectors,
and when the k vectors are in the xz plane.
<p><dt><code>run-yeven-zeven</code>, <code>run-yeven-zodd</code>, <code>run-yodd-zeven</code>, <code>run-yodd-zodd</code>, <code>run-te-yeven</code>, <code>run-te-yodd</code>, <code>run-tm-yeven</code>, <code>run-tm-yodd</code>
<dd>These <code>run</code>-like functions combine the
<code>yeven</code>/<code>yodd</code> constraints with
<code>zeven</code>/<code>zodd</code> or
<code>te</code>/<code>tm</code>. See also <code>run-parity</code>,
below.
<p><dt><code>(run-parity <i>p reset-fields band-func ...</i>)</code>
<dd>Like the <code>run</code> function, except that it takes two extra
parameters, a parity <code>p</code> and a boolean
(<code>true</code>/<code>false</code>) value
<code>reset-fields</code>. <code>p</code> specifies a parity
constraint, and should be one of the predefined variables:
<ul>
<li><code>NO-PARITY</code>: equivalent to <code>run</code>
<li><code>EVEN-Z</code> (or <code>TE</code>): equivalent to <code>run-zeven</code> or <code>run-te</code>
<li><code>ODD-Z</code> (or <code>TM</code>): equivalent to <code>run-zodd</code> or <code>run-tm</code>
<li><code>EVEN-Y</code> (like <code>EVEN-Z</code> but for y=0 plane)
<li><code>ODD-Y</code> (like <code>ODD-Z</code> but for y=0 plane)
</ul>
<p>It is possible to specify more than one symmetry constraint
simultaneously by adding them, e.g. <code>(+ EVEN-Z ODD-Y)</code>
requires the fields to be even through z=0 and odd through y=0. It is
an error to specify incompatible constraints (e.g. <code>(+ EVEN-Z
ODD-Z)</code>). <b>Important:</b> if you specify the z/y parity, the
dielectric structure (and the k vector) <b>must</b> be symmetric about the
z/y=0 plane, respectively.
<p>If <code>reset-fields</code> is <code>false</code>, the fields from
any previous calculation will be reused as the starting point from
this calculation, if possible; otherwise, the fields are reset to
random values. The ordinary <code>run</code> functions use a default
<code>reset-fields</code> of<code>true</code>. Alternatively,
<code>reset-fields</code> may be a string, the name of an HDF5 file to
load the initial fields from (as exported by
<code>save-eigenvectors</code>, <a href="#eigen-fields">below</a>).
<p><dt><code>(display-eigensolver-stats)</code>
<dd>Display some statistics on the eigensolver convergence; this function is
useful mainly for MPB developers in tuning the eigensolver.
</dl>
<p>Several band functions for outputting the eigenfields are defined
for your convenience, and are described in the <b>Band output
functions</b> section, below. You can also define your own band
functions, and for this purpose the functions described in the section
<b>Field manipulation functions</b>, below, are useful. A band
function takes the form:
<pre>
(define (<i>my-band-func</i> which-band)
<i>...do stuff here with band index which-band...</i>
)
</pre>
<p>Note that the output variable <code>freqs</code> may be used to
retrieve the frequency of the band (see above). Also, a global
variable <code>current-k</code> is defined holding the current k-point
vector from the <code>k-points</code> list.
<p>There are also some even lower-level functions that you can call,
although you should not need to do most of the time:
<dl>
<dt><code>(init-params <i>p reset-fields?</i>)</code>
<dd>Read the input variables and initialize the simulation in
preparation for computing the eigenvalues. The parameters are the
same as the first two parameters of <code>run-parity</code>.
This function <em>must</em> be called before any of the other
simulation functions below. (Note, however, that the <code>run</code>
functions all call <code>init-params</code>.)
<p><dt><code>(set-parity <i>p</i>)</code>
<dd>After calling <code>init-params</code>, you can change the
parity constraint without resetting the other parameters by
calling this function. Beware that this does not randomize the fields
(see below); you don't want to try to solve for, say, the TM
eigenstates when the fields are initialized to TE states from a
previous calculation.
<p><dt><code>(randomize-fields)</code>
<dd>Initialize the fields to random values.
<p><dt><code>(solve-kpoint <i>k</i>)</code>
<dd>Solve for the requested eigenstates at the Bloch wavevector
<code>k</code>.
</dl>
<h3><a name="find-k">The inverse problem: k as a function of frequency</a></h3>
<p>MPB's <code>(run)</code> function(s) and its underlying algorithms
compute the frequency <code>w</code> as a function of wavevector
<code>k</code>. Sometimes, however, it is desirable to solve the
inverse problem, for <code>k</code> at a given frequency
<code>w</code>. This is useful, for example, when studying coupling
in a waveguide between different bands at the same frequency
(frequency is conserved even when wavevector is not). One also uses
<code>k(w)</code> to construct wavevector diagrams, which aid in
understanding diffraction (e.g. negative-diffraction materials and
super-prisms). To solve such problems, therefore, we provide the
<code>find-k</code> function described below, which inverts
<code>w(k)</code> via a few iterations of Newton's method (using the
<a href="#group-vel">group velocity</a> <code>dw/dk</code>). Because
it employs a root-finding method, you need to specify bounds on
<code>k</code> and a <em>crude</em> initial guess (order of magnitude
is usually good enough).
<dl>
<dt><code>(find-k <i>p omega band-min band-max kdir
tol kmag-guess kmag-min kmag-max [band-func...]</i>)</code>
<dd>Find the wavevectors in the current geometry/structure for the
bands from <code><i>band-min</i></code> to
<code><i>band-max</i></code> at the frequency
<code><i>omega</i></code> along the <code><i>kdir</i></code> direction
in k-space. Returns a list of the wavevector magnitudes for each
band; the actual wavevectors are <code>(vector3-scale <i>magnitude</i>
(unit-vector3 <i>kdir</i>))</code>. The arguments of
<code>find-k</code> are:
<ul>
<li><code><i>p</i></code>: parity (same as first argument to
<code>run-parity</code>, <a href="#run">above</a>).
<li><code><i>omega</i></code>: the frequency at which to find the bands
<li><code><i>band-min</i></code>, <code><i>band-max</i></code>: the
range of bands to solve for the wavevectors of (inclusive).
<li><code><i>kdir</i></code>: the direction in k-space in which to
find the wavevectors. (The magnitude of <code><i>kdir</i></code> is
ignored.)
<li><code><i>tol</i></code>: the fractional tolerance with which to
solve for the wavevector; <code>1e-4</code> is usually sufficient.
(Like the <code>tolerance</code> input variable, this is only the
tolerance of the numerical iteration...it does not have anything to do
with e.g. the error from finite grid <code>resolution</code>.)
<li><code><i>kmag-guess</i></code>: an initial guess for the k
magnitude (along <code><i>kdir</i></code>) of the wavevector at
<code><i>omega</i></code>. Can either be a list (one guess for each
band from <code><i>band-min</i></code> to
<code><i>band-max</i></code>) or a single number (same guess for all
bands, which is usually sufficient).
<li><code><i>kmag-min</i></code>, <code><i>kmag-max</i></code>: a
range of k magnitudes to search; should be large enough to include the
correct k values for all bands.
<li><code><i>band-func</i></code>: zero or more <a
href="#band-funcs">band functions</a>, just as in <code>(run)</code>,
which are evaluated at the computed k points for each band.
</ul>
<p>The <code>find-k</code> routine also prints a line suitable for grepping:
<pre>
kvals: <i>omega</i>, <i>band-min</i>, <i>band-max</i>, <i>kdir-x</i>, <i>kdir-y</i>, <i>kdir-z</i>, <i>k magnitudes...</i>
</pre>
</dl>
<h3><a name="band-funcs">Band/output functions</a></h3>
<p>All of these are functions that, given a band index, output the
corresponding field or compute some function thereof (in the primitive
cell of the lattice). They are designed to be passed as band
functions to the <code>run</code> routines, although they can also be
called directly. See also the section on <a href="#field-norm">field
normalizations</a>.
<dl>
<dt><code>(output-hfield <i>which-band</i>)</code>
<dt><code>(output-hfield-x <i>which-band</i>)</code>
<dt><code>(output-hfield-y <i>which-band</i>)</code>
<dt><code>(output-hfield-z <i>which-band</i>)</code>
<dd>Output the magnetic (H) field for <code>which-band</code>; either
all or one of the components, respectively.
<p><dt><code>(output-dfield <i>which-band</i>)</code>
<dt><code>(output-dfield-x <i>which-band</i>)</code>
<dt><code>(output-dfield-y <i>which-band</i>)</code>
<dt><code>(output-dfield-z <i>which-band</i>)</code>
<dd>Output the electric displacement (D) field for
<code>which-band</code>; either all or one of the components,
respectively.
<p><dt><code>(output-efield <i>which-band</i>)</code>
<dt><code>(output-efield-x <i>which-band</i>)</code>
<dt><code>(output-efield-y <i>which-band</i>)</code>
<dt><code>(output-efield-z <i>which-band</i>)</code>
<dd>Output the electric (E) field for <code>which-band</code>; either
all or one of the components, respectively.
<p><dt><code>(output-hpwr <i>which-band</i>)</code>
<dd>Output the time-averaged magnetic-field energy density (hpwr =
|<b>H</b>|<sup><small>2</small></sup>) for <code>which-band</code>.
<p><dt><code>(output-dpwr <i>which-band</i>)</code>
<dd>Output the time-averaged electric-field energy density (dpwr =
epsilon*|<b>E</b>|<sup><small>2</small></sup>) for
<code>which-band</code>.
<p><dt><code>(fix-hfield-phase <i>which-band</i>)</code>
<dt><code>(fix-dfield-phase <i>which-band</i>)</code>
<dt><code>(fix-efield-phase <i>which-band</i>)</code>
<dd>Fix the phase of the given eigenstate in a canonical way based on
the given spatial field (see also <code>fix-field-phase</code>,
below). Otherwise, the phase is random; these functions also maximize
the real part of the given field so that one can hopefully just
visualize the real part. To fix the phase for output, pass one of
these functions to <code>run</code> before the corresponding output
function, e.g. <code>(run-tm fix-dfield-phase output-dfield-z)</code>
<p>Although we try to maximize the "real-ness" of the field, this has
a couple of limitations. First, the phase of the different field
components cannot, of course, be chosen independently, so an
individual field component may still be imaginary. Second, if you use
<code>mpbi</code> to take advantage of <a href="#inv-symmetry">inversion
symmetry</a> in your problem, the phase is mostly determined elsewhere
in the program; <code>fix-_field-phase</code> in that case only
determines the sign.
</dl>
<p>See also below for the <code>output-poynting</code> and
<code>output-tot-pwr</code> functions to output the Poynting vector
and the total electromagnetic energy density, respectively.
<p>Sometimes, you only want to output certain bands. For example,
here is a function that, given an band/output function like the ones
above, returns a new output function that only calls the first
function for bands with a large fraction of their energy in an object(s).
(This is useful for picking out defect states in supercell
calculations.)
<dl>
<dt><code>(output-dpwr-in-objects <i>band-func min-energy objects...</i>)</code>
<dd>Given a band function <code>band-func</code>, returns a new band
function that only calls <code>band-func</code> for bands having a
fraction of their electric-field energy greater than
<code>min-energy</code> inside the given objects (zero or more
geometric objects). Also, for each band, prints the fraction of their
energy in the objects in the following form (suitable for grepping):
<pre>
dpwr:, <i>band-index</i>, <i>frequency</i>, <i>energy-in-objects</i>
</pre>
</dl>
<p><code>output-dpwr-in-objects</code> only takes a single band
function as a parameter, but if you want it to call several band
functions, you can easily combine them into one with the following
routine:
<dl>
<dt><code>(combine-band-functions <i>band-funcs...</i>)</code>
<dd>Given zero or more band functions, returns a new band function
that calls all of them in sequence. (When passed zero parameters,