MPB光子晶體仿真軟件使用介紹

1、MPB 使用指南（部分）2006 年 12 月 01 日星期五MPB 使用指南 （部分 ）MPB 使用指南在這里， 我們會(huì)展示如何使用 MPB 進(jìn)行二維光子晶體能帶計(jì)算和輸出場(chǎng)分布圖的整個(gè)過(guò)程。你可以從中知道MPB 如何工作，也可以了解什么樣的東西可以用它來(lái) 實(shí)現(xiàn)。這里所列出的只是一部分，在MPB User Reference 里會(huì)有更詳細(xì)的內(nèi)容。在下一個(gè)專題， dataanalysis tutorial ，會(huì)有更多的例子，著重?cái)?shù)據(jù)的分析和可視化。ctl 文件在 MPB 的使用中， ctl 文件是不可缺少的，它的后綴是“ctl ，文件名類似” foo.ctl （你可以用你自己喜歡的名字代替fo

2、o ）。 ctl 文件包括了 所要研究的幾何結(jié)構(gòu)，要計(jì)算的本征矢量的數(shù)目，想要輸出的東西和其他你想要計(jì)算的東西。 ctl 是用腳本語(yǔ)言來(lái)寫的，所以它可以寫成一系列簡(jiǎn)單的命令，來(lái)設(shè)計(jì)幾何結(jié)構(gòu)等等。在這個(gè)文件中全部是用戶輸入，循環(huán)和其他必須的命令。不過(guò)不用擔(dān)心， 你不須要做一個(gè)真正的程序員， 因?yàn)檫@些語(yǔ)言都是比較簡(jiǎn)單的， 例如你可以不用按順序來(lái)輸入，不用理會(huì)空格，可以隨便插入說(shuō)明，也可以不理會(huì)其他默認(rèn)的設(shè)置。ctl 文件是執(zhí)行在libctl 庫(kù)上面的， 而 libctl 也是建立在 Scheme 語(yǔ)言上。 因此， 在一個(gè) ctl文件中有三種可能的命令和語(yǔ)法：1. Scheme 是由 MIT 開發(fā)出

3、來(lái)的一個(gè)強(qiáng)大的程序語(yǔ)言，它的語(yǔ)法很簡(jiǎn)單：所有的狀態(tài)量都是以下這個(gè)形式，（function arguments.）。我們要在GUN Guile 編譯器下運(yùn)行Scheme 。你不必學(xué)太多的 Scheme 來(lái)寫一個(gè)基本的 ctl 文件，你可以在需要的時(shí)候再去查找。當(dāng)然，有興趣的話，可以參考它 們的主頁(yè)。2. libctl 是我們用 Guile 編寫的一個(gè)庫(kù)，它是用來(lái)簡(jiǎn)化 Scheme 和科學(xué)計(jì)算軟件的接口。libctl 設(shè)置了一些基本的格式來(lái)實(shí)現(xiàn)用戶接口和定義大量有用的函數(shù)。具體參考其主頁(yè)。3. MPB 定義了全部的接口，用來(lái)實(shí)現(xiàn)光子能帶的計(jì)算。在manual 里，會(huì)著重說(shuō)明它的特點(diǎn)。如果你能去了解

4、一下libctl manual ， 你會(huì)獲益非淺， 特別是 libctl Basic User Experience那一節(jié)，這樣你就可以知道用戶接口是怎樣的， Scheme 語(yǔ)言大概是怎樣的(這個(gè)是很有用的)，還有一些有用的一般性質(zhì)。在這里我們就不再重復(fù)了。那就讓我們繼續(xù)。MPB 程序一般是用以下的命令來(lái)運(yùn)行：unix% mpb foo.ctl >& foo.out這樣，程序就會(huì)讀入 ctl 文件，并且執(zhí)行，保存數(shù)據(jù)在foo.out 這個(gè)文件里(在mpb-ctl /examples / 文件夾里有一些ctl 文件的例子)。當(dāng)然，你也可以直接輸入 mpb 命令，這樣就進(jìn)入了對(duì)話模式

5、，你可以繼續(xù)輸入命令然后直接看到結(jié)果。計(jì)算第一個(gè)能帶結(jié)構(gòu)我們第一個(gè)例子是計(jì)算由介質(zhì)棒構(gòu)成的四方晶格的二維能帶。在我們的 ctl 文件里，我們會(huì)首先指明要仿真的幾何結(jié)構(gòu)和參數(shù)，然后讓它運(yùn)行和輸出。所有的參數(shù)都有默認(rèn)值，每個(gè)參數(shù)對(duì)應(yīng)著一個(gè)Scheme 變量，所以我們只須指明哪些是需要修改的。(如果在guile 提示符下輸入命令： (help) ，程序會(huì)列出所有的參數(shù)變量和它們的信息)。其中一個(gè)參數(shù)是num-bands ，是指明在每個(gè)K 點(diǎn)要計(jì)算的能帶的數(shù)目。如果你在提示符下輸入 num-bands ，它會(huì)返回當(dāng)前值： 1 這個(gè)數(shù)值太小了，可以加大：(set! num-bands 8)這就是我們改變

6、參數(shù)的過(guò)程(如果你現(xiàn)在輸入是(與順序無(wú)關(guān))，設(shè)置我們想要計(jì)算能帶的它是一個(gè)三維矢量，默認(rèn)值是空的。我們把Gamma,X, M, 和 Gamma ：(set! k-points (list (vector3 0 0 0) ; Gamma(vector3 0.5 0 0) ; X(vector3 0.5 0.5 0) ; M(vector3 0 0 0) ; Gammanum-bands ，它會(huì)顯示8 )。下一步要做的K 點(diǎn)。這個(gè)參數(shù)由變量 k-points 來(lái)決定，K 點(diǎn)取在四方格子的簡(jiǎn)約布里淵區(qū)的轉(zhuǎn)角處，注意我們是怎樣建立一個(gè)list ，還有三個(gè)矢量vector3 ；我們可以分多行輸入，還可以

7、在分 號(hào)后加上注釋。一般上，我們也要計(jì)算三個(gè)方向上的 K 點(diǎn)的能帶，這樣才可以得到近似連續(xù)的能帶圖。我們可以調(diào)用一個(gè)函數(shù)來(lái)插入 K 點(diǎn)：(set! k-points (interpolate 4 k-points)這個(gè)函數(shù)可以在每個(gè)轉(zhuǎn)角間線性插入 4 個(gè) K 點(diǎn)；如果我們?cè)谔崾痉螺斎?k-points ，它會(huì)顯示以下的 16 個(gè)點(diǎn)：(#(0 0 0) #(0.1 0.0 0.0) #(0.2 0.0 0.0) #(0.3 0.0 0.0) #(0.4 0.0 0.0) #(0.5 00)#(0.5 0.1 0.0) #(0.5 0.2 0.0) #(0.5 0.3 0.0) #(0.5 0.

8、4 0.0) #(0.5 0.5 0)#(0.4 0.4 0.0) #(0.3 0.3 0.0) #(0.2 0.2 0.0) #(0.1 0.1 0.0) #(0 0 0)如上面描述的那樣，當(dāng) K 點(diǎn)(非歸一化的)定義在倒格矢上時(shí)，程序里全部的空間矢量都是在定義的格點(diǎn)基矢上， 而基矢是歸一在 basis-size 長(zhǎng)度的。 在 這種情況下，我們不必指明格點(diǎn)方向，因?yàn)槌绦蚰J(rèn)了格點(diǎn)的單位坐標(biāo)軸(例如，最常見(jiàn)的四方 / 立方晶格)(倒格點(diǎn)的矢量也是如此)。下面，我們將會(huì)看到如何 改變實(shí)空間格點(diǎn)基矢?，F(xiàn)在， 我們要設(shè)置晶體的幾何結(jié)構(gòu)， 所以我們先要指明在格點(diǎn)中心的元胞究竟包含什么東西，而變量 ge

9、ometry 就是我們需要設(shè)置的，它包括了幾個(gè)對(duì)象。在libctl 文件里提到，對(duì)象(object ) 是一個(gè)復(fù)雜的結(jié)構(gòu)數(shù)據(jù)類型， 它是由狀態(tài)量( (make type (property1 value1)(property2 value2) .) 產(chǎn)生的。 在幾何對(duì)象里， 包含了幾種變量： 圓柱棒 ( cylinders ) ,球體( spheres ) , 方塊( blocks ) , 橢球( ellipsoids )，在將來(lái)可能會(huì)有其他類型。現(xiàn)在，我們想要一個(gè)圓柱棒的四方晶格結(jié)構(gòu)，所以我們把一根介質(zhì)圓柱棒擺在格點(diǎn)中心：(set! geometry (list (make cylinder

10、(center 0 0 0) (radius 0.2) (height infinity)(material (make dielectric (epsilon 12)在這里， 我們?cè)O(shè)置了介質(zhì)棒的幾個(gè)性質(zhì)。 center 是在中心， radius 是 0.2 ， height 是無(wú)窮。另外一個(gè)性質(zhì)是material ，它自己也是一個(gè)對(duì)象 我們讓它的性質(zhì)是dielectric ，而且介電常數(shù) epsilon 是 12 。還有一個(gè)性質(zhì)我們可以設(shè)置，就是介質(zhì)棒的軸向，不過(guò)我們可以 保留默 認(rèn)，即是在 z 方向上。所有的幾何對(duì)象在表面上是三維的，但是我們要做的是二維仿真，唯一重要的是介質(zhì)棒和xy 平面

11、的交點(diǎn)。換句話說(shuō)，就是我們要設(shè)定晶體的維數(shù)。一般上，在我們定義計(jì)算單元的大小時(shí)， 這項(xiàng)工作已經(jīng)同時(shí)完成了。 計(jì)算單元的大小是由變量geometry-lattice 來(lái)決定的，它是 lattice 類的一個(gè)對(duì)象。我們可以設(shè)置一些維數(shù)為 no-size ，這樣就把這個(gè)維度去掉，減少了系統(tǒng)的維數(shù)：(set! geometry-lattice (make lattice (size 1 1 no-size)這里， 我們定義了一個(gè)1*1 的二維單元 (默認(rèn)是四方格子) 。 這個(gè)單元可以由變量resolution來(lái)離散化，默認(rèn)值是10 。這個(gè)數(shù)值對(duì)于二維的計(jì)算來(lái)說(shuō)比較小，所以我們可以加大：(set! re

12、solution 32)那么計(jì)算單元就會(huì)變成一個(gè)32*32 的計(jì)算網(wǎng)格。從計(jì)算的效率來(lái)說(shuō)，最好是把網(wǎng)格設(shè)成2的倍數(shù)，或者是一些小的數(shù)字的倍數(shù)(例如： 2 ， 3 ， 5 和 7 )。粗略估計(jì)，把resolution設(shè)成 8 以上才可以得到一個(gè)合理的精度。至此， 我們已經(jīng)完成了參數(shù)的設(shè)置， 其實(shí)還有幾個(gè)其他的參數(shù)， 但是在這里我們保留默認(rèn)值。我們現(xiàn)在就要計(jì)算能帶了。最簡(jiǎn)單的方法是輸入 (run) 。因?yàn)檫@是二維計(jì)算，所以我們可以分別計(jì)算 TM 模和 TE 模，這時(shí)，我們要輸入 (run-te) 和 (run-tm) 。運(yùn)行時(shí)，程序會(huì)顯示一大堆的輸出， 說(shuō)明程序的運(yùn)行狀態(tài)。這些輸出大多數(shù)是解釋，但

13、我們要特別指出的是其中的一個(gè)數(shù)據(jù)輸出：tefreqs: ， 13, 0.3, 0.3, 0, 0.424264, 0.372604, 0.540287, 0.644083, 0.81406,0.828135, 0.890673, 1.01328, 1.1124這些輸出可以允許你輕易地獲取數(shù)據(jù)并輸入畫圖軟件的數(shù)據(jù)表中。這些數(shù)據(jù)分別是K 點(diǎn)的序號(hào)， K 點(diǎn)的坐標(biāo)，能帶的頻率，它們由逗號(hào)分開。每一行都有前綴， tefreqs 代表 TE 能 帶 ,tmfreqs 代表 TM 能帶 ,freqs 代表一般的能帶(由 (run) 產(chǎn)生)。使用這些前綴，你可以通過(guò)程序grep 從輸出中抽取你想要的數(shù)據(jù)。例

14、如，如果你把輸出直接輸出到一個(gè)文件foo.out ，你可以在Unix 提示符 下輸入 grep tmfreqs foo.out 來(lái)抽取 TM 模的數(shù)據(jù)。注意輸出中包含一行解釋：tefreqs:, k index, kx, ky, kz, kmag/2pi, band 1, band 2, band 3, band 4, band 5,band 6, band 7, band 8這是每一列數(shù)據(jù)的名稱。另外一個(gè)運(yùn)行run 后的輸出是計(jì)算能帶中的帶隙。例如運(yùn)行(run-tm) 后的輸出包括了以下的帶隙輸出：Gap from band 1 (0.282623311147724) to band 2 (

15、0.419334798706834),38.9514660888911%Gap from band 4 (0.715673834754345) to band 5 (0.743682920649084),3.8385522650349%這些數(shù)據(jù)也存儲(chǔ)在一個(gè)變量gap-list ， 它包括了 (gap-percent ， gap-min ， gap-max) 。值得注意的是，在這些帶隙輸出中會(huì)存在著 “假帶隙 ” ，因?yàn)橛袃蓚€(gè)可能的原因：1. 如果兩個(gè)能帶交叉，一個(gè)假帶隙就會(huì)產(chǎn)生。因?yàn)槌绦蚴羌僭O(shè)能帶是不交叉的。這些假能帶一般都是很窄的 (< 1%) 。為了確認(rèn)起見(jiàn)，你必須考察模式的對(duì)稱性或

16、者計(jì)算交叉點(diǎn)附近的能帶。(即使交叉點(diǎn)在所計(jì)算的一個(gè)K 點(diǎn)上，也不會(huì)因?yàn)閿?shù)值的原因而產(chǎn)生一個(gè)準(zhǔn)確的簡(jiǎn) 并)2. 計(jì)算能帶時(shí)一般會(huì)取第一布里淵區(qū)邊界上的 K 點(diǎn)，但能帶也可能產(chǎn)生在布里淵區(qū)內(nèi)的 K點(diǎn)。為了確認(rèn)確實(shí)存在著帶隙(和大小)，你也要計(jì)算在布里淵區(qū)內(nèi)的 K 點(diǎn)的能帶。到這里，我們已經(jīng)可以計(jì)算能帶并抽取每個(gè)K 點(diǎn)的本征頻率了，但是我們還可以觀察場(chǎng)分布和介質(zhì)分布。不過(guò)我們需要輸出 HDF 文件。( HDF 是一個(gè)以二進(jìn)制存儲(chǔ)的多維的科學(xué)數(shù)據(jù)，可以用很多可視化軟件來(lái)顯示，我們可以用 HDF5 的格式輸出數(shù)據(jù)，它的后綴是« L ".h5 ）?！碑?dāng)你調(diào)用了其中一個(gè)run 函數(shù)是，計(jì)

17、算單元的介質(zhì)分布函數(shù)也會(huì)自動(dòng)輸出到文件“ epsilon.h5 。為了輸出場(chǎng)分布或者其他的信息，我們必須在 ”run 后面加上一些參數(shù)。例如：（run-tm output-efield-z）（run-te （output-at-kpoint （vector3 0.5 0 0） output-hfield-z output-dpwr）這樣就可以輸出所有K 點(diǎn)的 TM 模的電場(chǎng) （E） 的 z 分量，還有X 點(diǎn)的 TE 模的磁場(chǎng) （H） 的 z分量和能量強(qiáng)度（D）。輸出的文件名類似“e.k12.b03.z.te.h5: 表示TE卞H（.te）的第3個(gè)能帶（.b03）的第12個(gè)K點(diǎn)（.k12）的電場(chǎng)

18、（e）的z分量（.z）。每一個(gè)HDF5文件可以包含多維的數(shù)據(jù)，在這種情況下，它會(huì)包含場(chǎng)的實(shí)部和虛部（ z.r 和 z.i ），一般還包含了其他的數(shù)據(jù)（例如， output- efield 會(huì)在一個(gè)文件里輸出所有的分量）。參考Data AnalysisTutorial 。你也可以使用其他的幾個(gè)輸出函數(shù)，在User Reference 里列出，例如 output-dfield ，output-hpwr 和 output-dpwr-in-objects 。實(shí)際上，你可以隨便調(diào)用其他的函數(shù)，這樣你就可以做很多的事情， 而不局限于場(chǎng)分布的輸出， 你可以利用我們下面介紹的函數(shù)來(lái)實(shí)現(xiàn)能帶的分析。我們也可以直

19、接利用底層源代碼的功能來(lái)代替run 函數(shù)的功能，這樣就可以比較好地控制計(jì)算，例如在reference 一節(jié)里介紹的函數(shù)。關(guān)于單位原則上，在這個(gè)程序里，你可以定義任何你想要的選項(xiàng)。你可以輸入你自己的長(zhǎng)度和坐標(biāo)，例如 furlongs （英式單位），也可以有根據(jù)地設(shè)置格點(diǎn)基矢（在 geometry-lattice 里的 basis-size ）。然而，我們不提倡這樣做。Maxwell 方程組有一個(gè)重要的特點(diǎn)：尺度不變（與尺度無(wú)關(guān)）。如果你把尺寸擴(kuò)大到 10倍，方程的解也會(huì)擴(kuò)大到 10 倍（頻率會(huì)縮少到 1/10 ）。所以， 你可以一次性求解，然后如晶格常數(shù)，然把尺度換成你想要的單位?；谶@個(gè)原因，

20、我們經(jīng)常選擇一些基本的尺度，后以此作單位來(lái)計(jì)算其他的長(zhǎng)度。也就是說(shuō)，我們選擇了晶格常數(shù) a 作為單位，然后把所有長(zhǎng)度除以 a, 應(yīng)用到所有的物理量。 這個(gè)思想貫穿在前面和后面的例子， 也是 MPB 默認(rèn)的：晶格常數(shù)是一個(gè)單位，所有 的坐標(biāo)以此作參考。正如已經(jīng)提到的那樣， 程序里全部的三維矢量是定義在已經(jīng)規(guī)范化的格點(diǎn)基矢上。 即是每一個(gè)分量都是乘以相應(yīng)的基矢，然后合成相應(yīng)的矢量。所以在你設(shè)定的單位a 中，只有一般的坐標(biāo)尺寸，除此之外，所謂的長(zhǎng)度是沒(méi)有什么大的意義的。必須注意的是， 在上面提到的關(guān)于倒格點(diǎn)的方面有一個(gè)例外： 他們是定義在倒空間的， 所以如果定義是no-size ，那么他的倒格矢是2

21、*pi/a 。我們也設(shè)計(jì)了函數(shù)用來(lái)轉(zhuǎn)換不同庫(kù)的矢量。程序解出來(lái)的本征頻率是以c/a 作為單位的， c 是真空光速， a 是晶格常數(shù)(所以，相應(yīng)的真空波長(zhǎng)就是a 除以本征頻率)。三角晶格的能帶作為第二個(gè)例子， 我們會(huì)計(jì)算一個(gè)由介質(zhì)棒構(gòu)成的三角晶格的 TM 能帶。 我們只要通過(guò)變量geometry-lattice 來(lái)改變格點(diǎn)就行了。我們要把兩個(gè)基矢設(shè)成與x 軸夾角都為 30 度，且分別在 x 軸的上面和下面，這樣就改變了默認(rèn)的情況。(set! geometry-lattice (make lattice (size 1 1 no-size)(basis1 (/ (sqrt 3) 2) 0.5)(b

22、asis2 (/ (sqrt 3) 2) -0.5)basic3 不用改變， 保留它作為 z 軸。 注意到 Scheme 為我們提供了一般的運(yùn)算符號(hào)和函數(shù)，但也保留了自己的風(fēng)格： 使用前綴。 在這里， 我們只是定義了矢量的方向， 并沒(méi)有定義長(zhǎng)度，依然也是默認(rèn)的單位長(zhǎng)度，這是無(wú)傷大雅的。因?yàn)槿蔷Ц竦暮?jiǎn)約布里淵區(qū)不同于四方晶格，所以我們要相應(yīng)改變要計(jì)算的 k( k-points ) :(set! k-points (list (vector3 0 0 0) ; Gamma(vector3 0 0.5 0) ; M(vector3 (/ -3) (/ 3) 0) ; K(vector3 0 0 0

23、) ; Gamma (set! k-points (interpolate 4 k-points)注意，這些矢量定義了新的倒格矢，所以 k 點(diǎn)坐標(biāo)也是在倒格矢坐標(biāo)上的，并非是在正交的坐標(biāo)軸上。(在Scheme 中， (/ 3) 也相當(dāng)于 (/ 1 3) 或者 1/3 )。所有其他參數(shù)( geometry,num-bands, 和 grid-size )可以保留與前面的相同，這樣我們就可以馬上運(yùn)行(run-tm) 來(lái)計(jì)算能帶。結(jié)果出來(lái)后，我們可以看見(jiàn) TM 模的帶隙比前面的四方晶格要寬：Gap from band 1 (0.275065617068082) to band 2 (0.446289

24、918847647), 47.4729292989213%注: 這是美國(guó)麻省理工學(xué)院的 MPB 的 tutorial 的一部分 , 由本人翻譯,如有出入 ,請(qǐng)登陸其主頁(yè)以獲取正確的信息。b.mpb(GPL) 提供能帶計(jì)算的平面波展開法，包含麥克斯韋方程組的變換展開。詳細(xì)見(jiàn)源代碼。c.h5utils(GPL) 是一個(gè)可視化的組件，提供生成二維圖和三維可視化文件， .v5d 和 .vtk ，可以分別使用vis5d 和 paraview 軟件，具體這些程序哪里下載和安裝，請(qǐng)自行搜索或者 Email 我。d.qtiplot(GPL) ， Qt 界面的繪圖軟件，功能齊全。e.libctl(GPL) 是一

25、個(gè)關(guān)于常用晶體格子的函數(shù)庫(kù)，提供晶格和格點(diǎn)的 “描述” ，向上支持調(diào)用，典型的軟件有 MPB 和 MEEP 。f.fftw 是廣泛使用的離散傅立葉變換庫(kù)之一，類似有GSL ， LAPACK 和私有的 Intel MKL 等。要注意的是，以上軟件使用 GNU GPL 自由軟件協(xié)議發(fā)行。好的，現(xiàn)在默認(rèn)你已經(jīng)安裝了以上軟件。如果不能很好安裝或者想使用最新版本， 請(qǐng)到 MPB 官方網(wǎng)站 下載對(duì)應(yīng)軟件的最新版本， 按照安裝提示自行解壓 安裝。現(xiàn)在要運(yùn)行測(cè)試了，假設(shè)目錄位置為 mpb/ ，則編寫以下文件：main.ctl 和 run.sh123456789101112131415161718 192021

26、 222324252627282930 31 3233 343536 37;main.#!/bin/basctlh mpb- (set!main.ctl >num-bands 8);&這里設(shè)置八個(gè)main.out h帶，是每一個(gè)k5topng -S 3點(diǎn)計(jì)算帶的個(gè)epsilon.h5數(shù) (set!h5tov5dgeometry-lattepsilon.h5ice (makeh5tovtklattice (size12345678epsilon.h51 1 no-size)9101112131mpb-data(basis1 (/4151617-r -m 3 -n(sqrt 3) 2)320.5)epsilon.h5h5ls epsilon.h5h5topng epsilon.h5 :data-new h5tov5d epsilon.h5 :data-new grep Gap(basis2 (/(sqrt 3) 2)-0.5) ; 這里設(shè)置晶格的組成結(jié)構(gòu)， 晶格基矢為 (1 1no-size) ， 因?yàn)槭瞧矫娴脑?;這里是元包的m