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
|
; Test suite for MPB. This file runs MPB for a variety of cases,
; and compares it against known results from previous versions. If the
; answers aren't sufficiently close, it exits with an error.
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Some general setup and utility routines first:
(set! tolerance 1e-9) ; use a low tolerance to get consistent results
; keep track of some error statistics:
(define min-err infinity)
(define max-err 0)
(define max-abs-err 0)
(define sum-err 0)
(define sum-abs-err 0)
(define num-err 0)
; function to check if two results are sufficently close:
(define-param check-tolerance 1e-3)
(define (almost-equal? x y)
(if (> (abs x) 1e-3)
(let ((err (/ (abs (- x y)) (* 0.5 (+ (abs x) (abs y)))))
(abserr (abs (- x y))))
(set! min-err (min min-err err))
(set! max-err (max max-err err))
(set! max-abs-err (max max-abs-err abserr))
(set! num-err (+ num-err 1))
(set! sum-err (+ sum-err err))
(set! sum-abs-err (+ sum-abs-err abserr))))
(or
(< (abs (- x y)) (* 0.5 check-tolerance (+ (abs x) (abs y))))
(and (< (abs x) 1e-3) (< (abs (- x y)) 1e-3))))
; Convert a list l into a list of indices '(1 2 ...) of the same length.
(define (indices l)
(if (null? l)
'()
(cons 1 (map (lambda (x) (+ x 1)) (indices (cdr l))))))
; Check whether the freqs returned by a run (all-freqs) match correct-freqs.
(define (check-freqs correct-freqs)
(define (check-freqs-aux fc-list f-list ik)
(define (check-freqs-aux2 fc f ib)
(if (not (almost-equal? fc f))
(error "check-freqs: k-point " ik " band " ib " is "
f " instead of " fc)))
(if (= (length fc-list) (length f-list))
(map check-freqs-aux2 fc-list f-list (indices f-list))
(error "check-freqs: wrong number of bands at k-point " ik)))
(if (= (length correct-freqs) (length all-freqs))
(begin
(map check-freqs-aux correct-freqs all-freqs (indices all-freqs))
(print "check-freqs: PASSED\n"))
(error "check-freqs: wrong number of k-points")))
; checks whether list X and list Y are almost equal
(define (check-almost-equal X Y)
(if (fold-left (lambda (x y) (and x y)) true (map almost-equal? X Y))
(print "check-almost-equal: PASSED\n")
(error "check-almost-equal: FAILED\n" X Y)))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(if (not (using-mpi?)) ; MPI code currently doesn't support 1d systems
(begin
; Use a lower tolerance for the 1d cases, since it is cheap; otherwise,
; the Bragg-sine case perennially causes problems.
(set! tolerance (/ tolerance 10000))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; First test: a simple 1d Bragg mirror:
(print
"**************************************************************************\n"
" Test case: 1d quarter-wave stack.\n"
"**************************************************************************\n"
)
(set! geometry (list (make cylinder (material (make dielectric (epsilon 9.0)))
(center 0) (axis 1)
(radius infinity) (height 0.25))))
(set! k-points (interpolate 4 (list (vector3 0 0 0) (vector3 0.5 0 0))))
(set! grid-size (vector3 32 1 1))
(set! num-bands 8)
(let ((correct-freqs '((0.0 0.666384282528928 0.666667518031881 1.33099337665679 1.33336087672839 1.99161980173015 2.00024302642565 2.6450937852083) (0.0574931097997446 0.608787973096119 0.724373977055668 1.2735966394103 1.39097797892597 1.93584941706027 2.05633589084836 2.59377005482918) (0.113352271592983 0.552761237115776 0.780756082813777 1.21672824612875 1.44856381682192 1.87756171149864 2.11568713317692 2.53392260905135) (0.164802158007456 0.501201811770545 0.832987385247272 1.16417795261225 1.50250265212664 1.82314371802974 2.17224997337128 2.47690271740002) (0.205536502065922 0.460405353660882 0.8747810869492 1.12220337812548 1.54664162749105 1.77873062187023 2.22033216569935 2.42854622021095) (0.22245099319197 0.443470718153721 0.89236781440248 1.10455805248897 1.56579373692671 1.75948932731251 2.24248043100853 2.40631190043537))))
(run-tm)
(check-freqs correct-freqs)
(run-te)
(check-freqs correct-freqs))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Second test: a less-simple 1d Bragg mirror, consisting of a sinusoidally
; varying dielectric index (see also bragg-sine.ctl):
(print
"**************************************************************************\n"
" Test case: 1d sinusoidal Bragg mirrors.\n"
"**************************************************************************\n"
)
(let ((pi (* 4 (atan 1)))) ; 3.14159...
(define (eps-func p)
(make dielectric (index (+ 2 (cos (* 2 pi (vector3-x p)))))))
(set! default-material (make material-function (material-func eps-func))))
(set! k-points (interpolate 9 (list (vector3 0 0 0) (vector3 0.5 0 0))))
(set! grid-size (vector3 32 1 1))
(set! num-bands 8)
(run-tm)
(check-freqs '((0.0 0.460648275218079 0.542427739817 0.968586586361068 1.01617062660004 1.48337333794755 1.48386583676098 1.96763768634053) (0.0231424888788602 0.454293827620926 0.548894464606965 0.958360293771355 1.02641972845204 1.45913774058533 1.50811173706602 1.94948827608003) (0.0462090787082577 0.439084320039512 0.564452582533386 0.938267239403995 1.04658249131212 1.43467215953867 1.53260886786687 1.92547645623973) (0.0691102159040715 0.420015536860369 0.584144708540882 0.915749746086578 1.06922078451581 1.41022612312792 1.55710983305974 1.90110410737745) (0.0917238673443258 0.399480768741237 0.60565316278248 0.892483769139075 1.09266692571585 1.38580829656309 1.58161025389159 1.87664635327897) (0.113863114251085 0.378567678109852 0.628026307077094 0.868958449423545 1.11644695897109 1.3614327391061 1.60610415581645 1.85216493367767) (0.135212561098235 0.357979431767488 0.650804940461757 0.845399327507094 1.14036437789844 1.33712597279413 1.63058104368944 1.82768870377694) (0.155193837618337 0.338479349371984 0.673671347997618 0.822003372220962 1.16428444411062 1.31294698795159 1.65501798584125 1.80324819761608) (0.172679547014293 0.321289854992633 0.696193190997784 0.799137496016584 1.18800188305626 1.28906327514362 1.67934881274624 1.77891315965595) (0.185502873728775 0.308627942742938 0.717034417258616 0.778101847588796 1.21076249234602 1.26620930567935 1.70326114132094 1.75499775017206) (0.19041596916807 0.303766163770654 0.728590051056375 0.766483258741006 1.22530241888082 1.25163924613769 1.72099865592632 1.73725912578794)))
(set! default-material air) ; don't screw up later tests
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(set! tolerance (* tolerance 10000))
)
) ; if (not (using-mpi?))
; Test get-dominant-planwave
(print
"**************************************************************************\n"
" Test case: get-dominant-planewave.\n"
"**************************************************************************\n"
)
(set! geometry '())
(set! k-points (list (vector3 0.4 0 0)))
(set! geometry-lattice (make lattice (size 1 0.1 no-size)))
(set! num-bands 8)
(set! resolution 100)
(run-te)
(check-almost-equal (list (vector-ref (get-dominant-planewave 1) 0)
(vector-ref (get-dominant-planewave 2) 0)
(vector-ref (get-dominant-planewave 3) 0))
(list 0.4 -0.6 1.4))
;; Restore initial values
(set! resolution 32)
(set! geometry-lattice (make lattice (size 1 1 no-size)))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Square lattice of dielectric rods in air.
(print
"**************************************************************************\n"
" Test case: Square lattice of dielectric rods in air.\n"
"**************************************************************************\n"
)
(set! geometry (list
(make cylinder (material (make dielectric (epsilon 11.56)))
(center 0 0) (radius 0.2) (height infinity))))
(set! k-points (interpolate 4 (list (vector3 0) (vector3 0.5)
(vector3 0.5 0.5 0) (vector3 0))))
(set! grid-size (vector3 32 32 1))
(set! num-bands 8)
(run-te)
(check-freqs '((0.0 0.561945334557385 0.780842812331228 0.780846098412592 0.924371683163949 1.00803985129792 1.00804088731802 1.09858623956652) (0.0897997798994076 0.560209896469707 0.767785350870472 0.782438231618754 0.912603067546361 1.00888633849104 1.00949114803277 1.12115474664388) (0.178852773321546 0.553300983639649 0.732522222968759 0.786672449175983 0.890463105855642 1.01385729138557 1.02078013833592 1.11300646463179) (0.266123594110195 0.534865350840295 0.689376138533469 0.792049548816113 0.872925911023147 1.02090314736827 1.04285404891473 1.10753383632769) (0.349588020537016 0.494779669995165 0.658508223213229 0.796523676734899 0.862683642813515 1.02898182530566 1.07043694197557 1.10003079693163) (0.413345586803412 0.44462291351553 0.648672217588454 0.798265471046059 0.859327463282679 1.03311041395003 1.09580214986717 1.09741680977047) (0.424298362126911 0.44898279373346 0.644548562548272 0.802551148059984 0.854853159661748 0.98983998591444 1.0583076027211 1.11679157683255) (0.455353836640927 0.461160511700901 0.633230253238328 0.814782424365419 0.834250744952672 0.934557961696952 1.01274334726807 1.12367050454956) (0.478427219924468 0.501912134249683 0.617449702512683 0.784489069111248 0.833710248078536 0.906540901997911 0.967733581444411 1.1263583529554) (0.495604528258234 0.556867480138193 0.601739184739058 0.720584183505631 0.858431813314025 0.897771643330533 0.92661062844362 1.12761853567986) (0.503605801531701 0.594344024613055 0.59435663953178 0.679388741318615 0.883664739567964 0.895781768469706 0.895786874556172 1.12801276576574) (0.474670773656218 0.549917137140298 0.608031881017845 0.745034531026827 0.848354751485403 0.895854164182374 0.945890643267912 1.12708247876502) (0.373253932819571 0.543225052446258 0.646250803261438 0.817389372081735 0.830457382271798 0.896516227567732 1.01760084433849 1.12197791661875) (0.252369333718973 0.551073128162047 0.700804337123788 0.797551602017551 0.899585769590366 0.903389586014203 1.08791421563576 1.10099000357684) (0.126940001741566 0.558853684776988 0.755651182474766 0.785047154907932 0.909884122992341 0.968812657472461 1.04813919691317 1.11332213355591) (0.0 0.561945334407731 0.780842812401179 0.780846098384021 0.924371683195103 1.00803985129514 1.00804088729554 1.09858622997793)))
(run-tm)
(check-freqs '((0.0 0.550336075492761 0.561337783494192 0.561339793996441 0.822948013585295 0.868841613389014 0.965325380929893 1.08937760109445) (0.0651416429431535 0.525003695654183 0.561884878000649 0.586306893982008 0.823535977020037 0.86734647097958 0.954632345309176 1.05687677175817) (0.127664746818747 0.493649813919317 0.563322581824471 0.617311826620811 0.822736862854134 0.863529964061771 0.92430556209052 1.03882168368046) (0.184046514295327 0.461592151644409 0.565121689802102 0.651396635807459 0.810838592407438 0.85898411247184 0.892949156581082 1.03589199660828) (0.227778015582534 0.433360712061513 0.566595731512458 0.689044651877333 0.778791875677485 0.855425103295414 0.879398204894174 1.0387432995892) (0.245808974576747 0.420657338406186 0.56716328782128 0.720091820469093 0.747202991063479 0.854090458576806 0.877011871859037 1.04079703189466) (0.249298773583017 0.427307623035431 0.560220797848635 0.718025037316212 0.756122437610169 0.855013330357386 0.877106455296421 1.02916965498279) (0.258693333199135 0.445540869703809 0.543385384840764 0.711697647242147 0.779350400764203 0.858135812263486 0.877357018542926 1.00092092539592) (0.270962613899996 0.470322466744217 0.524052873443864 0.701396122540261 0.810963308664218 0.864393838558208 0.877683037781189 0.965108093819696) (0.281613643862846 0.493207050642694 0.508823220586706 0.68996057277849 0.846784814196112 0.87459516255628 0.878068231243395 0.926580909129451) (0.285905779127161 0.502981364580489 0.502983097838737 0.684476386658726 0.874359380527121 0.883317372585053 0.883317410406254 0.892993349560143) (0.276088805107277 0.491352389673393 0.508683441988827 0.692579899252015 0.839723261134077 0.85643142702439 0.907218621287156 0.907347258076394) (0.240239034783233 0.479230170582085 0.523498859086802 0.685361670368096 0.829265933706741 0.840450969771508 0.910791381218572 0.941595252480836) (0.175240325825588 0.488533364028353 0.541560691215046 0.647509833455847 0.83001959727719 0.850443167942857 0.922657201444286 0.983924775906362) (0.0915258938978901 0.516393376876295 0.555923804221748 0.60121063231939 0.824561228580069 0.865290375816656 0.948411582724858 1.03526242687651) (0.0 0.550336075497221 0.561337783491839 0.561339793993807 0.822948013590914 0.868841613389918 0.965325380904055 1.08937790223254)))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(print
"****************************************************************************\n"
" Test case: square lattice of magneto-electric rods in air.\n"
"****************************************************************************\n"
)
; Test case provided by Ling Lu, who checked the results against COMSOL for the
; hermitian-eps case.
(set! num-bands 10)
(set! geometry-lattice (make lattice (size 1 1 no-size)))
(set! resolution 32)
(set! k-points
(interpolate 2 (list (vector3 0 0) (vector3 0.5 0) (vector3 0.5 0.5) (vector3 0 0))))
(set! default-material air)
(set! geometry
(list (make cylinder (center 0 0 0) (radius 0.11) (height infinity)
(material (make medium-anisotropic
(epsilon-diag 15 15 15) (epsilon-offdiag 0 0 0)
(mu-diag 14 14 1) (mu-offdiag (if (has-hermitian-eps?)
0+12.4i 1.1) 0 0))))))
(run-tm)
(check-freqs (if (has-hermitian-eps?) '((0.0 0.463294018556493 0.578605647982777 0.651041639469486 0.89507841938234 0.906562954329028 0.971861717371308 1.00315281180164 1.01424590331597 1.02482341106969) (0.129380123044635 0.461952809408746 0.581559682819551 0.650968421346652 0.851508658392909 0.899132837287078 0.948338264572998 0.986837587655266 1.02012949457141 1.11135500583651) (0.24061857781097 0.455581797026516 0.594368290992035 0.650900996467599 0.74835731456284 0.900869187217702 0.951667861016453 0.985036605138704 1.0516822808838 1.12149556246633) (0.29226915952455 0.447982173597163 0.613801580856078 0.650943532698264 0.687698304986519 0.901519360659148 0.957067156979669 0.984771165606936 1.08285810218551 1.11061141488713) (0.300271180585828 0.469567031887115 0.597375716765173 0.674305570222573 0.692027497826549 0.897950544588226 0.926840934886075 0.969149696783797 0.998326180698363 1.10650926174014) (0.316662188469676 0.51370370227899 0.588292652694574 0.689922957143738 0.730249942797529 0.845838260405215 0.859587865047183 0.913738876925751 0.986600858452482 1.13738513559847) (0.325124793328178 0.529201233604094 0.601833025098704 0.702804083485205 0.779409771583802 0.781916228087728 0.801633572250502 0.910192402752782 0.985981745195422 1.15351027935509) (0.299048781303725 0.498217631142972 0.581862947853374 0.670216351165454 0.79230544712122 0.801008335876475 0.901133548767124 0.934783461321372 0.987891734100305 1.1274295455455) (0.179169782471429 0.464237623484746 0.581862665371791 0.652256768573312 0.857854382579286 0.865493306970474 0.914866356757855 0.983643482419635 1.03320702306621 1.12806207145723) (0.0 0.463294040286175 0.578605723955974 0.651041643962349 0.895078379737811 0.906562937848727 0.97186159963616 1.00315283598839 1.01424612908797 1.02482337953567)) '((0.0 0.26588586999867 0.356853038485526 0.368990587212674 0.503897363414621 0.506551995618629 0.539911334424117 0.635682859764425 0.645852865039137 0.66007402686939) (0.123330481479188 0.278293838663916 0.357366805753587 0.369461558917164 0.503797813078635 0.506641848145437 0.538218354695777 0.635573939560424 0.645599854632284 0.659512681920588) (0.194658392275999 0.331458695546422 0.361679789289526 0.377947926726132 0.502693772691889 0.506755517926571 0.527230780909839 0.633336678239672 0.644216769995567 0.653387932908076) (0.208804178427075 0.348977997743993 0.365278518063831 0.443856727877108 0.477902469744969 0.506844674794979 0.509075958831114 0.617529235852054 0.639721519490117 0.646946211798357) (0.210676648364611 0.351106119979908 0.364916521569197 0.453356051624939 0.500502931100297 0.507018635289265 0.51012548275551 0.616826993693061 0.642495569786502 0.646905514765383) (0.214081523585675 0.353217483148627 0.365002039172578 0.471976957273301 0.507165103362517 0.51079005950879 0.563895397008847 0.615053559332354 0.646182218771027 0.648364548586666) (0.215649021943222 0.353602936491205 0.365269847254897 0.480748866729077 0.507335500482168 0.512534454027105 0.607921694331747 0.619786381474157 0.648912790283471 0.650328008406889) (0.210879518955412 0.354630556068409 0.361246363774902 0.439488506919655 0.507022083817614 0.508252106018903 0.543820254199556 0.620599301809527 0.646605526208855 0.651154928068886) (0.162353699928896 0.295876053296023 0.356214325457118 0.372161020104797 0.50425817165436 0.506735257227901 0.537356240859274 0.633231503464085 0.645934609764122 0.65717004035603) (0.0 0.265892439896466 0.356853509421417 0.368991043876895 0.50389707509854 0.506552257076724 0.539913534063763 0.635682631384305 0.645872010464673 0.660017515310061))))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Using the targeted solver to find a defect state in a 5x5 triangular
; lattice of rods.
(print
"**************************************************************************\n"
" Test case: 3x3 triangular lattice of rods in air, dipole defect states.\n"
"**************************************************************************\n"
)
(if (not force-mu?) ; targeted solver doesn't handle mu yet
(begin
(set! geometry-lattice (make lattice (size 3 3 1)
(basis1 (/ (sqrt 3) 2) 0.5)
(basis2 (/ (sqrt 3) 2) -0.5)))
(set! k-points (list (vector3 0 0.5 0))) ; K
(set! geometry (list
(make cylinder (material (make dielectric (epsilon 12)))
(center 0 0) (radius 0.2) (height infinity))))
(set! geometry (geometric-objects-lattice-duplicates geometry))
(set! geometry (append geometry
(list (make cylinder (center 0 0 0)
(radius 0.33) (height infinity)
(material (make dielectric (epsilon 12)))))))
(set! grid-size (vector3 (* 16 5) (* 16 5) 1))
(set! num-bands 2)
(set! target-freq 0.35)
(run-tm)
(let ((ct-save check-tolerance))
(set! check-tolerance (* ct-save 10))
(check-freqs '((0.33627039929402 0.338821383601027)))
(set! check-tolerance ct-save))
(set! target-freq 0)))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(print
"**************************************************************************\n"
" Test case: fcc lattice of air spheres in dielectric.\n"
"**************************************************************************\n"
)
(set! geometry-lattice (make lattice
(basis1 0 1 1)
(basis2 1 0 1)
(basis3 1 1 0)))
(set! k-points (interpolate 1 (list
(vector3 0 0.5 0.5) ; X
(vector3 0 0.625 0.375) ; U
(vector3 0 0.5 0) ; L
(vector3 0 0 0) ; Gamma
(vector3 0 0.5 0.5) ; X
(vector3 0.25 0.75 0.5) ; W
(vector3 0.375 0.75 0.375)))) ; K
(set! geometry (list (make sphere (center 0) (radius 0.5) (material air))))
(set! default-material (make dielectric (epsilon 11.56)))
(set! grid-size (vector3 16 16 16))
(set! mesh-size 5)
(set! num-bands 10)
(run)
(check-freqs '((0.3703065073388731 0.3720847951594916 0.3825201413648734 0.383600351309727 0.4939893048065466 0.5123606721164143 0.5223073978045091 0.523821730596807 0.5951030481294071 0.6640811083906616) (0.36835908725434985 0.37748681319649185 0.3849917491239541 0.3874383238948237 0.4718945404212059 0.5065458897880816 0.5237249784688186 0.531091150473746 0.6102926625627835 0.6497145516818431) (0.3575328801327589 0.380949914017912 0.39278926172097595 0.4014539592513005 0.43934988081812454 0.49515765158554836 0.5274783178063797 0.5408918418169267 0.6355155083750688 0.6446208129570911) (0.32374150485195863 0.3316528408482444 0.39820951394232795 0.40120559265452554 0.4628679442744884 0.5142922894295575 0.5332819309482675 0.5461160730637494 0.6310458540642596 0.6442359011147937) (0.30752940384296407 0.3088914413473654 0.387976766798534 0.3899131086118616 0.4916056724921297 0.5363288685004943 0.5366126847030026 0.538984533419493 0.6246809592147906 0.6291240037206337) (0.1797554271422153 0.1804214571518261 0.47340167673207667 0.4759683186484289 0.5034596764371356 0.5357700886824948 0.5381441770965257 0.5403446090795936 0.6240734738388297 0.627676429794933) (0.0 0.0 0.5182834534581704 0.5211243503626508 0.5211324776053197 0.5439419138927141 0.5439517918334954 0.5465047798955615 0.6128700951828456 0.6139704268282182) (0.2072575103849531 0.20784064700695526 0.47322381806952135 0.47526776006597526 0.5068725664421999 0.5268739560935114 0.5303663197703882 0.5320228430098722 0.6046890975929905 0.6539550359503241) (0.3703065491828492 0.37208484798757957 0.3825202709398482 0.3836004905486825 0.4939893802033828 0.5123607939536433 0.5223076744430113 0.5238218663879185 0.5951031328487622 0.6640854671485138) (0.37230930969457554 0.3771541648897106 0.38596761425016163 0.392561609652401 0.46356851069458777 0.5020256500434837 0.505399132586034 0.549256352020148 0.6215017804667959 0.6392031720356646) (0.3737093998668942 0.38640413833892384 0.38649906634404985 0.41020828777713325 0.4352292457022003 0.4908086507266613 0.49251835949389555 0.5665833567717682 0.627970994624438 0.6540781805440787) (0.36426053680366377 0.3827541756966786 0.39001099138782824 0.4059668434549568 0.43825785151278884 0.49360218014587165 0.5057362637756401 0.5592192482696754 0.6348958863356928 0.6450932995888251) (0.3591221748572636 0.3802582312189367 0.3912758255354526 0.4026805821435226 0.4380053809311872 0.49516016649470224 0.5291218538334365 0.5416787182268407 0.637583231663347 0.6456239632717847)))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(print
"**************************************************************************\n"
" Test case: simple cubic lattice with anisotropic dielectric.\n"
"**************************************************************************\n"
)
(set! geometry-lattice (make lattice))
(set! default-material air)
(set! k-points (list (vector3 0) (vector3 0.5)
(vector3 0.5 0.5) (vector3 0.5 0.5 0.5)))
(set! grid-size (vector3 16 16 16))
(set! mesh-size 5)
(define hi-all (make dielectric (epsilon 12)))
(define hi-x (make dielectric-anisotropic (epsilon-diag 12 1 1)))
(define hi-y (make dielectric-anisotropic (epsilon-diag 1 12 1)))
(define hi-z (make dielectric-anisotropic (epsilon-diag 1 1 12)))
(set! geometry
(list (make block (center 0) (size 0.313 0.313 1) (material hi-z))
(make block (center 0) (size 0.313 1 0.313) (material hi-y))
(make block (center 0) (size 1 0.313 0.313) (material hi-x))
(make block (center 0) (size 0.313 0.313 0.313)
(material hi-all))))
(set! num-bands 3)
(run)
(check-freqs '((0.0 0.0 0.546634963647193) (0.259951207183097 0.259951258670471 0.444658075510584) (0.300692330235002 0.345673935773948 0.497692646240215) (0.362782432577119 0.362782474830613 0.502236936705387)))
(print
"*******************************************************************************\n"
" Test case: group velocity in simple cubic lattice with anisotropic dielectric.\n"
"*******************************************************************************\n"
)
(set! k-points (list (vector3 0.12 0.34 0.41)))
(run)
(let ((v (compute-group-velocity-component (vector3 1 2 3))))
(check-almost-equal v (map (lambda (b) (compute-1-group-velocity-component (vector3 1 2 3) b))
(arith-sequence 1 1 num-bands)))
(check-almost-equal v '(0.202671224992983 0.310447990695762 -0.0480795046912859)))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(if (not (using-mpi?)) ; currently do not support transformed-overlap with MPI
(begin
(print
"**************************************************************************\n"
" Test case: symmetry transformed overlaps & inversion/mirror eigenvalues.\n"
"**************************************************************************\n"
)
(set! resolution 16)
(set! num-bands 6)
(set! default-material air)
; define a simple geometry with inversion and z-mirror
(set! geometry-lattice
(make lattice (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 0.25) (material (make dielectric (epsilon 13) )))))
; compare inversion eigenvalues at k = [1,1,1]/2 against precomputed values
(set! k-points (list (vector3 0.5 0.5 0.5))) ; little group includes inversion
(let ((Winv (matrix3x3 (vector3 -1 0 0) (vector3 0 -1 0) (vector3 0 0 -1)))
(Wmz (matrix3x3 (vector3 1 0 0) (vector3 0 1 0) (vector3 0 0 -1)))
(w (vector3 0 0 0)) )
(run)
; check inversion eigenvalues (inversion as operation = {Winv|w})
(let ((symeigs-inv (compute-symmetries Winv w)))
(check-almost-equal (map real-part symeigs-inv) '(+1 +1 +1 -1 -1 -1))
(check-almost-equal (map imag-part symeigs-inv) '( 0 0 0 0 0 0))
)
; compare with run-zeven and run-zodd at k = 0 [must be 0 due to https://github.com/NanoComp/mpb/issues/114]
(set! k-points (list (vector3 0 0 0)))
(run-zeven) ; even z-parity check
(let ((symeigs-zeven (compute-symmetries Wmz w)))
(check-almost-equal (map real-part symeigs-zeven) '(+1 +1 +1 +1 +1 +1))
(check-almost-equal (map imag-part symeigs-zeven) '( 0 0 0 0 0 0))
)
(run-zodd) ; odd z-parity check
(let ((symeigs-zodd (compute-symmetries Wmz w)))
(check-almost-equal (map real-part symeigs-zodd) '(-1 -1 -1 -1 -1 -1))
(check-almost-equal (map imag-part symeigs-zodd) '( 0 0 0 0 0 0))
)
) ; let ((Winv ...
) ; begin
) ; if (not (using-mpi?))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(display-eigensolver-stats)
(print "Relative error ranged from " min-err " to " max-err
", with a mean of " (/ sum-err num-err) "\n")
(print "Absolute error ranged to " max-abs-err
", with a mean of " (/ sum-abs-err num-err) "\n")
(print "PASSED all tests.\n")
|
