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

MPB使用指南（部分）2006年12月01日星期五MPB使用指南（部分）MPB使用指南在这里，我们会展示如何使用MPB进行二维光子晶体能带计算和输出场分布图的整个过程你可以从中知道MPB如何工作，也可以了解什么样的东西可以用它来实现这里所列出的只是一部分，在MPBUserReference里会有更详细的内容在下一个专题，dataanalysistutorial，会有更多的例子，着重数据的分析和可视化ctl文件在MPB的使用中，ctl文件是不可缺少的，它的后缀是“ctl，”文件名类似foo.ctl（你可以用你自己喜欢的名字代替foo）ctl文件包括了所要研究的几何结构，要计算的本征矢量的数目，想要输出的东西和其他你想要计算的东西ctl是用脚本语言来写的，所以它可以写成一系列简单的命令，来设计几何结构等等在这个文件中全部是用户输入，循环和其他必须的命令不过不用担心，你不须要做一个真正的程序员，因为这些语言都是比较简单的，例如你可以不用按顺序来输入，不用理会空格，可以随便插入说明，也可以不理会其他默认的设置ctl文件是执行在libctl库上面的，而libctl也是建立在Scheme语言上因此，在一个ctl文件中有三种可能的命令和语法：Scheme是由MIT开发出来的一个强大的程序语言，它的语法很简单：所有的状态量都是以下这个形式，（functionarguments...）。

我们要在GUNGuile编译器下运行Scheme你不必学太多的Scheme来写一个基本的ctl文件，你可以在需要的时候再去查找当然，有兴趣的话，可以参考它们的主页1. libctl是我们用Guile编写的一个库，它是用来简化Scheme和科学计算软件的接口libctl设置了一些基本的格式来实现用户接口和定义大量有用的函数具体参考其主页2. MPB定义了全部的接口，用来实现光子能带的计算在manual里，会着重说明它的特点那一节，这样你就可以知道用户接口是怎样的，Scheme语言大概是怎样的(这个是很有用的)，还有一些有用的一般性质在这里我们就不再重复了那就让我们继续MPB程序一般是用以下的命令来运行：unix%mpbfoo.ctl>&foo.out这样，程序就会读入ctl文件，并且执行，保存数据在foo.out这个文件里(在mpb-ctl/examples/文件夹里有一些ctl文件的例子)当然，你也可以直接输入mpb命令，这样就进入了对话模式，你可以继续输入命令然后直接看到结果计算第一个能带结构我们第一个例子是计算由介质棒构成的四方晶格的二维能带在我们的ctl文件里，我们会首先指明要仿真的几何结构和参数，然后让它运行和输出。

所有的参数都有默认值，每个参数对应着一个Scheme变量，所以我们只须指明哪些是需要修改的如果在guile提示符下输入命令：(help)，程序会列出所有的参数变量和它们的信息)其中一个参数是num-bands，是指明在每个K点要计算的能带的数目如果你在提示符下输入num-bands，它会返回当前值：1——这个数值太小了，可以加大：(set!num-bands8)这就是我们改变参数的过程(如果你现在输入是(与顺序无关)，设置我们想要计算能带的它是一个三维矢量，默认值是空的我们把Gamma,X,M,和Gamma：(set!k-points(list(vector3000);Gamma(vector30.500);Xnum-bands，它会显示8)下一步要做的K点这个参数由变量k-points来决定，K点取在四方格子的简约布里渊区的转角处，(vector30.50.50);M(vector3000)));Gamma注意我们是怎样建立一个list，还有三个矢量vector3；我们可以分多行输入，还可以在分号后加上注释一般上，我们也要计算三个方向上的K点的能带，这样才可以得到近似连续的能带图我们可以调用一个函数来插入K点：(set!k-points(interpolate4k-points))这个函数可以在每个转角间线性插入4个K点；如果我们在提示符下输入k-points，它会显示以下的16个点：(#(000)#(0.10.00.0)#(0.20.00.0)#(0.30.00.0)#(0.40.00.0)#(0.500)#(0.50.10.0)#(0.50.20.0)#(0.50.30.0)#(0.50.40.0)#(0.50.50)#(0.40.40.0)#(0.30.30.0)#(0.20.20.0)#(0.10.10.0)#(000))如上面描述的那样，当K点(非归一化的)定义在倒格矢上时，程序里全部的空间矢量都是在定义的格点基矢上，而基矢是归一在basis-size长度的。

在这种情况下，我们不必指明格点方向，因为程序默认了格点的单位坐标轴(例如，最常见的四方/立方晶格)(倒格点的矢量也是如此)下面，我们将会看到如何改变实空间格点基矢现在，我们要设置晶体的几何结构，所以我们先要指明在格点中心的元胞究竟包含什么东西，而变量geometry就是我们需要设置的，它包括了几个对象在libctl文件里提到，对象(object)是一个复杂的结构数据类型，它是由状态量((maketype(property1value1)(property2value2)...))产生的在几何对象里，包含了几种变量：圆柱棒(cylinders),球体(spheres),方块(blocks),椭球(ellipsoids)，在将来可能会有其他类型现在，我们想要一个圆柱棒的四方晶格结构，所以我们把一根介质圆柱棒摆在格点中心：(set!geometry(list(makecylinder(center000)(radius0.2)(heightinfinity)(material(makedielectric(epsilon12))))))在这里，我们设置了介质棒的几个性质center是在中心，radius是0.2，height是无穷。

另外一个性质是material，它自己也是一个对象——我们让它的性质是dielectric，而且介电常数epsilon是12还有一个性质我们可以设置，就是介质棒的轴向，不过我们可以保留默认，即是在z方向上所有的几何对象在表面上是三维的，但是我们要做的是二维仿真，唯一重要的是介质棒和xy平面的交点换句话说，就是我们要设定晶体的维数一般上，在我们定义计算单元的大小时，这项工作已经同时完成了计算单元的大小是由变量geometry-lattice来决定的，它是lattice类的一个对象我们可以设置一些维数为no-size，这样就把这个维度去掉，减少了系统的维数：(set!geometry-lattice(makelattice(size11no-size)))这里，我们定义了一个1*1的二维单元(默认是四方格子)这个单元可以由变量resolution来离散化，默认值是10这个数值对于二维的计算来说比较小，所以我们可以加大：(set!resolution32)那么计算单元就会变成一个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提示符下输入greptmfreqsfoo.out来抽取TM模的数据注意输出中包含一行解释：tefreqs:,kindex,kx,ky,kz,kmag/2pi,band1,band2,band3,band4,band5,band6,band7,band8这是每一列数据的名称。

另外一个运行run后的输出是计算能带中的带隙例如运行(run-tm)后的输出包括了以下的带隙输出：Gapfromband1(0.282623311147724)tobandGapfrombandtoband5(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-tmoutput-efield-z)(run-te(output-at-kpoint(vector30.500)output-hfield-zoutput-dpwr))这样就可以输出所有K点的TM模的电场(E)的z分量，还有X点的TE模的磁场(H)的z分量和能量强度(D)输出的文件名类似“”，表示TE模(.te)的第3个能带(.b03)的第12个K点(.k12)的电场(e)的z分量(.z)每一个HDF5文件可以包含多维的数据，在这种情况下，它会包含场的实部和虚部(z.r和z.i)，一般还包含了其他的数据(例如，output-efield会在一个文件里输出所有的分量)参考DataAnalysisTutorial你也可以使用其他的几个输出函数，在UserReference里列出，例如output-dfield，output-hpwr和output-dpwr-in-objects实际上，你可以随便调用其他的函数，这样你就可以做很多的事情，而不局限于场分布的输出，你可以利用我们下面介绍的函数来实现能带的分析。

我们也可以直接利用底层源代码的功能来代替run函数的功能，这样就可以比较好地控制计算，例如在reference一节里介绍的函数关于单位原则上，在这个程序里，你可以定义任何你想要的选项你可以输入你自己的长度和坐标，例如furlongs(英式单位)，也可以有根据地设置格点基矢(在geometry-lattice里的basis。