MPB 光子晶体仿真软件使用介绍

MPB使用指南(部分)
2006年12月01日星期五
MPB使用指南(部分)
MPB使用指南
在这里，我们会展示如何使用MPB进行二维光子晶体能带计算和输出场分布图的整个过程。

你可以从中知道MPB如何工作，也可以了解什么样的东西可以用它来实现。

这里所列出的只是一部分，在MPB User Reference里会有更详细的内容。

在下一个专题，data analysis tutorial，会有更多的例子，着重数据的分析和可视化。

ctl文件
在MPB的使用中，ctl文件是不可缺少的，它的后缀是“ctl”，文件名类似foo.ctl（你可以用你自己喜欢的名字代替foo）。

ctl文件包括了所要研究的几何结构，要计算的本征矢量的数目，想要输出的东西和其他你想要计算的东西。

ctl是用脚本语言来写的，所以它可以写成一系列简单的命令，来设计几何结构等等。

在这个文件中全部是用户输入，循环和其他必须的命令。

不过不用担心，你不须要做一个真正的程序员，因为这些语言都是比较简单的，例如你可以不用按顺序来输入，不用理会空格，可以随便插入说明，也可以不理会其他默认的设置。

ctl文件是执行在libctl库上面的，而libctl也是建立在Scheme语言上。

因此，在一个ctl 文件中有三种可能的命令和语法：
1. Scheme 是由MIT开发出来的一个强大的程序语言，它的语法很简单：所有的状态量都是以下这个形式，(function arguments...) 。

我们要在GUN Guile编译器下运行Scheme。

你不必学太多的Scheme来写一个基本的ctl文件，你可以在需要的时候再去查找。

当然，有兴趣的话，可以参考它们的主页。

2. libctl 是我们用Guile编写的一个库，它是用来简化Scheme和科学计算软件的接口。

libctl设置了一些基本的格式来实现用户接口和定义大量有用的函数。

具体参考其主页。

3. MPB 定义了全部的接口，用来实现光子能带的计算。

在manual里，会着重说明它的特点。

如果你能去了解一下libctl manual，你会获益非浅，特别是libctl Basic User Experience 那一节，这样你就可以知道用户接口是怎样的，Scheme语言大概是怎样的（这个是很有用的），还有一些有用的一般性质。

在这里我们就不再重复了。

那就让我们继续。

MPB程序一般是用以下的命令来运行：
unix% mpb foo.ctl >& foo.out
这样，程序就会读入ctl文件，并且执行，保存数据在foo.out这个文件里（在mpb-ctl / examples / 文件夹里有一些ctl文件的例子）。

当然，你也可以直接输入mpb命令，这样就进入了对话模式，你可以继续输入命令然后直接看到结果。

计算第一个能带结构
我们第一个例子是计算由介质棒构成的四方晶格的二维能带。

在我们的ctl文件里，我们会首先指明要仿真的几何结构和参数，然后让它运行和输出。

所有的参数都有默认值，每个参数对应着一个Scheme变量，所以我们只须指明哪些是需要修改的。

（如果在guile提示符下输入命令：(help) ，程序会列出所有的参数变量和它们的信息）。

其中一个参数是num-bands ，是指明在每个K点要计算的能带的数目。

如果你在提示符下输入num-bands，它会返回当前值：1 ——这个数值太小了，可以加大：(set! num-bands 8)
这就是我们改变参数的过程（如果你现在输入num-bands，它会显示8）。

下一步要做的是（与顺序无关），设置我们想要计算能带的K点。

这个参数由变量k-points来决定，它是一个三维矢量，默认值是空的。

我们把K点取在四方格子的简约布里渊区的转角处，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))) ; Gamma
注意我们是怎样建立一个list，还有三个矢量vector3；我们可以分多行输入，还可以在分号后加上注释。

一般上，我们也要计算三个方向上的K点的能带，这样才可以得到近似连续的能带图。

我们可以调用一个函数来插入K点：
(set! k-points (interpolate 4 k-points))
这个函数可以在每个转角间线性插入4个K点；如果我们在提示符下输入k-points，它会显示以下的16个点：
(#(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 0 0)
#(0.5 0.1 0.0) #(0.5 0.2 0.0) #(0.5 0.3 0.0) #(0.5 0.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))
如上面描述的那样，当K点（非归一化的）定义在倒格矢上时，程序里全部的空间矢量都是在定义的格点基矢上，而基矢是归一在basis-size长度的。

在这种情况下，我们不必指明格点方向，因为程序默认了格点的单位坐标轴（例如，最常见的四方/立方晶格）（倒格点的矢量也是如此）。

下面，我们将会看到如何改变实空间格点基矢。

现在，我们要设置晶体的几何结构，所以我们先要指明在格点中心的元胞究竟包含什么东西，而变量geometry就是我们需要设置的，它包括了几个对象。

在libctl文件里提到，对象（object）是一个复杂的结构数据类型，它是由状态量（(make type (property1 value1) (property2 value2) ...)）产生的。

在几何对象里，包含了几种变量：圆柱棒（cylinders）, 球体（spheres）, 方块（blocks）,椭球（ellipsoids），在将来可能会有其他类型。

现在，我们想要一个圆柱棒的四方晶格结构，所以我们把一根介质圆柱棒摆在格点中心：
(set! geometry (list (make cylinder
(center 0 0 0) (radius 0.2) (height infinity)
(material (make dielectric (epsilon 12))))))
在这里，我们设置了介质棒的几个性质。

center是在中心，radius是0.2，height是无穷。

另外一个性质是material，它自己也是一个对象——我们让它的性质是dielectric，而且
介电常数epsilon是12。

还有一个性质我们可以设置，就是介质棒的轴向，不过我们可以保留默认，即是在z方向上。

所有的几何对象在表面上是三维的，但是我们要做的是二维仿真，唯一重要的是介质棒和xy平面的交点。

换句话说，就是我们要设定晶体的维数。

一般上，在我们定义计算单元的大小时，这项工作已经同时完成了。

计算单元的大小是由变量geometry-lattice来决定的，它是lattice类的一个对象。

我们可以设置一些维数为no-size，这样就把这个维度去掉，减少了系统的维数：
(set! geometry-lattice (make lattice (size 1 1 no-size)))
这里，我们定义了一个1*1的二维单元（默认是四方格子）。

这个单元可以由变量resolution 来离散化，默认值是10。

这个数值对于二维的计算来说比较小，所以我们可以加大：(set! resolution 32)
那么计算单元就会变成一个32*32的计算网格。

从计算的效率来说，最好是把网格设成2的倍数，或者是一些小的数字的倍数（例如：2，3，5和7）。

粗略估计，把resolution 设成8以上才可以得到一个合理的精度。

至此，我们已经完成了参数的设置，其实还有几个其他的参数，但是在这里我们保留默认值。

我们现在就要计算能带了。

最简单的方法是输入(run)。

因为这是二维计算，所以我们可以分别计算TM模和TE模，这时，我们要输入(run-te)和(run-tm)。

运行时，程序会显示一大堆的输出，说明程序的运行状态。

这些输出大多数是解释，但我们要特别指出的是其中的一个数据输出：
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
这些输出可以允许你轻易地获取数据并输入画图软件的数据表中。

这些数据分别是K点的序号，K点的坐标，能带的频率，它们由逗号分开。

每一行都有前缀，tefreqs代表TE能带,tmfreqs代表TM能带,freqs代表一般的能带（由(run) 产生）。

使用这些前缀，你
可以通过程序grep从输出中抽取你想要的数据。

例如，如果你把输出直接输出到一个文件foo.out，你可以在Unix提示符下输入grep tmfreqs foo.out来抽取TM模的数据。

注意输出中包含一行解释：
tefreqs:, k index, kx, ky, kz, kmag/2pi, band 1, band 2, band 3, band 4, band 5, band 6, band 7, band 8
这是每一列数据的名称。

另外一个运行run后的输出是计算能带中的带隙。

例如运行(run-tm)后的输出包括了以下的带隙输出：
Gap from band 1 (0.282623311147724) to band 2 (0.419334798706834), 38.9514660888911%
Gap from band 4 (0.715673834754345) to band 5 (0.743682920649084), 3.8385522650349%
这些数据也存储在一个变量gap-list，它包括了(gap-percent，gap-min，gap-max)。

值得注意的是，在这些带隙输出中会存在着“假带隙”，因为有两个可能的原因：
1. 如果两个能带交叉，一个假带隙就会产生。

因为程序是假设能带是不交叉的。

这些假能带一般都是很窄的(< 1%)。

为了确认起见，你必须考察模式的对称性或者计算交叉点附近的能带。

（即使交叉点在所计算的一个K点上，也不会因为数值的原因而产生一个准确的简并）
2. 计算能带时一般会取第一布里渊区边界上的K点，但能带也可能产生在布里渊区内的K 点。

为了确认确实存在着带隙（和大小），你也要计算在布里渊区内的K点的能带。

到这里，我们已经可以计算能带并抽取每个K点的本征频率了，但是我们还可以观察场分布和介质分布。

不过我们需要输出HDF文件。

（HDF是一个以二进制存储的多维的科学数据，可以用很多可视化软件来显示，我们可以用HDF5的格式输出数据，它的后缀是“.h5”）。

当你调用了其中一个run函数是，计算单元的介质分布函数也会自动输出到文件“epsilon.h5”。

为了输出场分布或者其他的信息，我们必须在run后面加上一些参数。

例如：
(run-tm output-efield-z)
(run-te (output-at-kpoint (vector3 0.5 0 0) output-hfield-z output-dpwr))
这样就可以输出所有K点的TM模的电场(E)的z分量，还有X点的TE模的磁场(H)的z 分量和能量强度(D)。

输出的文件名类似“e.k12.b03.z.te.h5”，表示TE模(.te)的第3个能带(.b03)的第12个K点(.k12)的电场(e)的z分量(.z)。

每一个HDF5文件可以包含多维的数据，在这种情况下，它会包含场的实部和虚部（z.r和z.i），一般还包含了其他的数据（例如，output- efield 会在一个文件里输出所有的分量）。

参考Data Analysis Tutorial。

你也可以使用其他的几个输出函数，在User Reference里列出，例如output-dfield，output-hpwr和output-dpwr-in-objects。

实际上，你可以随便调用其他的函数，这样你就可以做很多的事情，而不局限于场分布的输出，你可以利用我们下面介绍的函数来实现能带的分析。

我们也可以直接利用底层源代码的功能来代替run函数的功能，这样就可以比较好地控制计算，例如在reference 一节里介绍的函数。

关于单位
原则上，在这个程序里，你可以定义任何你想要的选项。

你可以输入你自己的长度和坐标，例如furlongs（英式单位），也可以有根据地设置格点基矢（在geometry-lattice里的basis-size）。

然而，我们不提倡这样做。

Maxwell方程组有一个重要的特点：尺度不变（与尺度无关）。

如果你把尺寸扩大到10倍，方程的解也会扩大到10倍（频率会缩少到1/10）。

所以，你可以一次性求解，然后把尺度换成你想要的单位。

基于这个原因，我们经常选择一些基本的尺度，如晶格常数，然
后以此作单位来计算其他的长度。

也就是说，我们选择了晶格常数a作为单位，然后把所有长度除以a,应用到所有的物理量。

这个思想贯穿在前面和后面的例子，也是MPB默认的：晶格常数是一个单位，所有的坐标以此作参考。

正如已经提到的那样，程序里全部的三维矢量是定义在已经规范化的格点基矢上。

即是每一个分量都是乘以相应的基矢，然后合成相应的矢量。

所以在你设定的单位a中，只有一般的坐标尺寸，除此之外，所谓的长度是没有什么大的意义的。

必须注意的是，在上面提到的关于倒格点的方面有一个例外：他们是定义在倒空间的，所以如果定义是no-size，那么他的倒格矢是2*pi/a。

我们也设计了函数用来转换不同库的矢量。

程序解出来的本征频率是以c/a作为单位的，c是真空光速，a是晶格常数（所以，相应的真空波长就是a除以本征频率）。

三角晶格的能带
作为第二个例子，我们会计算一个由介质棒构成的三角晶格的TM能带。

我们只要通过变量geometry-lattice来改变格点就行了。

我们要把两个基矢设成与x轴夹角都为30度，且分别在x轴的上面和下面，这样就改变了默认的情况。

(set! geometry-lattice (make lattice (size 1 1 no-size)
(basis1 (/ (sqrt 3) 2) 0.5)
(basis2 (/ (sqrt 3) 2) -0.5)))
basic3不用改变，保留它作为z轴。

注意到Scheme为我们提供了一般的运算符号和函数，但也保留了自己的风格：使用前缀。

在这里，我们只是定义了矢量的方向，并没有定义长度，依然也是默认的单位长度，这是无伤大雅的。

因为三角晶格的简约布里渊区不同于四方晶格，所以我们要相应改变要计算的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))) ; Gamma
(set! k-points (interpolate 4 k-points))
注意，这些矢量定义了新的倒格矢，所以k点坐标也是在倒格矢坐标上的，并非是在正交的坐标轴上。

（在Scheme中，(/ 3)也相当于(/ 1 3) 或者1/3）。

所有其他参数（geometry,num-bands,和grid-size）可以保留与前面的相同，这样我们就可以马上运行(run-tm)来计算能带。

结果出来后，我们可以看见TM模的带隙比前面的四方晶格要宽：
Gap from band 1 (0.275065617068082) to band 2 (0.446289918847647), 47.4729292989213%
注: 这是美国麻省理工学院的MPB的tutorial的一部分,由本人翻译,如有出入,请登陆其主页以获取正确的信息。

b.mpb(GPL)提供能带计算的平面波展开法，包含麦克斯韦方程组的变换展开。

详细见源代码。

c.h5utils(GPL)是一个可视化的组件，提供生成二维图和三维可视化文件，.v5d和.vtk，可以分别使用vis5d和paraview软件，具体这些程序哪里下载和安装，请自行搜索或者Email我。

d.qtiplot(GPL)，Qt界面的绘图软件，功能齐全。

e.libctl(GPL)是一个关于常用晶体格子的函数库，提供晶格和格点的“描述”，向上支持调用，典型的软件有MPB和MEEP。

f.fftw是广泛使用的离散傅立叶变换库之一，类似有GSL，LAPACK和私有的Intel MKL等。

要注意的是，以上软件使用GNU GPL自由软件协议发行。

好的，现在默认你已经安装了以上软件。

如果不能很好安装或者想使用最新版本，请到MPB官方网站下载对应软件的最新版本，按照安装提示自行解压安装。

现在要运行测试了，假设目录位置为mpb/，则编写以下文件：
main.ctl 和run.sh
123456789101112131415161 7181920212223242526272829 3031323334353637 ;-------main.
ctl----------
---(set!
num-bands 8);
这里设置八个
带，是每一个k
点计算带的个
数(set!
geometry-latt
ice (make
lattice (size
1 1 no-size)
(basis1 (/
(sqrt 3) 2)
0.5)
(basis2 (/
(sqrt 3) 2)
-0.5)));这里
设置晶格的组
成结构，晶格基
矢为(1 1
no-size)，因为
是平面的原因;
这里是元包的
12345678
9101112131
4151617
#!/bin/bas
h mpb
main.ctl >
&
main.out h
5topng -S 3
epsilon.h5
h5tov5d
epsilon.h5
h5tovtk
epsilon.h5
mpb-data
-r -m 3 -n
32
epsilon.h5
h5ls
epsilon.h5
h5topng
epsilon.h5
:data-new
h5tov5d
epsilon.h5
:data-new
grep Gap
两个基矢，注意语法(set! geometry (list (make
cylinder (center 0 0 0) (radius 0.2) (height infinity) (material (make dielectric (epsilon 12))))));设置格点位置放置
的晶系，有方格，球，椭圆，圆柱，锥体，台体;这里设置圆柱，中心(0 0 0)，半径0.2，高为无穷大，绝缘介质中，介电常数e=12;具
体的设置见官
方文档(set! k-points (list (vector3 0 0 0) ; Gamma
(vector3 0 0.5 0) ; M (vector3 (/
-3) (/ 3) 0) ; K
(vector3 0 0 0))) ; Gamma;设置布里渊区倒空间
的k点模样，由于是三角格子，倒空间也是三
角格子，;注意对称性的关系main.out g rep tmfreqs main.out > main.tm.da t grep tefreqs main.out > main.te.da t
可以减少范围
(set!
k-points
(interpolate 4
k-points));设
置k点，使每个
k点间隔4个点
内差值细分
(set!
resolution
32);设置分辨
率，也就是精确
度(run-tm
(output-at-kp
oint (vector3
(/ -3) (/ 3) 0)
fix-efield-ph
ase
output-efield
-z));计算TM
模式在矢量方
向的数值
(run-te);计
算TE模式
在/mpb目录下输入命令：
$ chmod 755 run.sh //该为可执行文件
$ ./run.sh //执行
等待结果，那么就可以得到图片和数据。

至于怎样绘制main.tm.dat和main.te.dat的数据？首先删除这些文件的第一行，替换删除“tmfreqs:,”，除去“、”替换为空格。

好了，打开qtiplot把数据读取进去，划出来就可以了。

这里粗糙画画。