The MPB Manual
From AbInitio
Contents
|
Introduction
| |
| MIT Photonic Bands |
| Download |
| Release notes |
| MPB manual |
| Introduction |
| Installation |
| User Tutorial |
| Data Analysis Tutorial |
| User Reference |
| Developer Information |
| Acknowledgements |
| License and Copyright |
MIT Photonic Bands is a software package to compute definite-frequency eigenstates of Maxwell's equations in periodic dielectric structures. Its primary intended application is the study of photonic crystals: periodic dielectric structures exhibiting a band gap in their optical modes, prohibiting propagation of light in that frequency range. However, it could also be easily applied to compute other optical dispersion relations and eigenstates (e.g. for conventional waveguides such as fiber-optic cables).
This manual assumes that the reader is familiar with concepts from solid-state physics such as eigenstates, band structures, and Bloch's theorem. We also do not attempt (much) to instruct the reader on photonic crystals or other optical applications for which this code might be useful. For an excellent introduction to all of these topics in the context of photonic crystals, see the book Photonic Crystals: Molding the Flow of Light, by J. D. Joannopoulos, R. D. Meade, and J. N. Winn (Princeton, 1995).
Some of the main design goals we were thinking about when we developed this package were the following (see also the feature list at the MPB home page):
- Fully vectorial, three-dimensional calculations for arbitrary Bloch wavevectors. (The only approximation is the spatial discretization, or equivalently the planewave cutoff.)
- Flexible interface. Readable, extensible, scriptable...see also the libctl design goals.
- Parallel. (Can run on a single-processor machine, but is also supports parallel machines with MPI.)
- "Targeted" eigensolver: find modes nearest to a specified frequency, not just the lowest-frequency bands. (For defect calculations.)
- Leverage existing software (LAPACK, BLAS, FFTW, HDF, MPI, GUILE...).
- Modularity. The eigensolver, Maxwell's equations, user interface, and so on, should be oblivious to each other as much as possible. This way, they can be debugged separately, combined in various ways, replaced, used in other programs...all the usual benefits of modular design.
- Take advantage of inversion symmetry in the dielectric function, but don't require it. (This means that we have to handle both real and complex fields.)
Frequency-Domain vs. Time-Domain
There are two common computational approaches to studying dielectric structures such as photonic crystals: frequency-domain and time-domain. We feel that each has its own place in a researcher's toolbox, and each has unique advantages and disadvantages (for a more detailed discussion, see our online textbook, appendix D). MPB is frequency-domain; that is, it does a direct computation of the eigenstates and eigenvalues of Maxwell's equations (using a planewave basis). Each field computed has a definite frequency. In contrast, time-domain techniques (as in our Meep software) iterate Maxwell's equations in time; the computed fields have a definite time (at each time step) but not a definite frequency per se. It seems worthwhile to say a few words about each method, and to explain why we want a frequency-domain code. (Our group has both time-domain and frequency-domain software, and we use both techniques extensively.)
Time-domain methods are well-suited to computing things that involve evolution of the fields, such as transmission and resonance decay-time calculations. They can also be used to calculate band structures and for finding resonant modes, by looking for peaks in the Fourier transform of the system response to some input. The main advantage of this is that you get all the frequencies (peaks) at once from a calculation involving propagation of a single field. There are several disadvantages to this technique, however. First, it is hard to be confident that you have found all of the states—you may have coupled weakly to some state by accident, or two states may be close in frequency and appear as a single peak; this is especially problematic in higher-order resonant-cavity and waveguide calculations. Second, in the Fourier transform, the frequency resolution is inversely related to the simulation time; to get 10 times the resolution you must run your simulation 10 times as long (although matters are improved by using more sophisticated signal-processing methods such as Harminv's). Third, the time-step size must be proportional to the spatial-grid size for numerical stability; thus, if you double the spatial resolution, you must double the number of time steps (the length of your simulation), even if you are looking at states with the same frequency as before. Fourth, you only get the frequencies of the states; to get the eigenstates themselves (so that you can see what the modes look like and do calculations with them), you must run the simulation again, once for each state that you want, and for a time inversely proportional to the frequency-spacing between adjacent states (i.e. a long time for closely-spaced states).
In contrast, frequency-domain methods like those in MPB are in many ways better-suited to calculating band structures and eigenstates. (Here, we consider the case of an iterative eigensolver like the one in MPB, which iteratively improves approximate eigenstates. Dense solvers, which factorize the matrix directly, are impractical for large problems because of the huge size of the matrix, and because they compute many more eigenvectors than are desired. In iterative methods, the operator is only applied to individual vectors and is never itself computed explicitly.) First, you don't have to worry about missing states—even closely-spaced modes will appear as two eigenvalues in the result. Second, the error in the frequency in an iterative eigensolver typically decays exponentially with the number of iterations, so the number of iterations is logarithmic in the desired tolerance. Third, the number of iterations typically remains almost constant even as you increase the resolution (the work for each iteration increases, of course, but that happens in time-domain too). Fourth, you get both the frequencies and the eigenstates at the same time, so you can look at the modes immediately (even closely-spaced ones).
A traditional disadvantage of frequency-domain methods was that you had to compute all of the lowest eigenstates, up to the desired one, even if you didn't care about the lower ones. This was especially problematic in defect calculations, in which a large supercell is used, because in that case the lower bands are "folded" many times in the Brillouin zone. Thus, you often had to compute a large number of bands in order to get to the one you wanted (incurring large costs both in time and in storage). These disadvantages disappear to some degree in MPB, however, with the advent of its "targeted" eigensolver—with it, you can solve directly for the localized defect states (i.e the states in the band gap) without computing the lower bands. However, the targeted eigensolver method used in MPB still has poor convergence, so time-domain methods such as Meep still have an advantage here.
History
In order to shed some light on the past development and design decisions of the Photonic-Bands package, I thought I'd write a few paragraphs about its history. Read on to enjoy my narcissistic ramblings.
Many different people have written codes to compute the modes of periodic dielectric materials (although we don't know of any others that have been freely released). I have had experience with several programs written within our group, and this experience has guided the design of MIT Photonic-Bands.
The first program of this sort that I came in contact with, and the code that was used in our group until the development of this package, was initially written around 1990 (in Fortran 77) by R. D. Meade. As this software grew organically over time, several problems became apparent. First, the input format was inflexible (difficult to add new features without breaking old simulations), sensitive to whitespace and other formatting, and required repeated entry of information that was often the same from file to file. Often, pre- and post-processing steps were required using additional scripts and tools. Second, the many parts of the program had become intertwined, lacking modularity or a clear flow of control; this made it difficult to follow or modify substantially. Parallelizing it, or removing constraints that it imposed like inversion symmetry, or even replacing the input format seemed impractical. Besides, even reading code with variables named gxgzco (in common block cabgv) (honest!) is a mind-altering experience.
After an initial experience in the Spring of 1996 at writing a code based on a wavelet, rather than a planewave, basis (which turned out not to be practical), I set out to write a replacement for our Fortran eigensolver. My main aims at this time (Fall 1996) were a more flexible and powerful input format and a code that would be amenable to parallelization. I succeeded in achieving a working code with similar convergence to the old code (after some pain), but I discovered several things. I had lots of fun learning to use lex and yacc (Yet Another Compiler-Compiler) to make a flexible, C-like input format with variables and other advanced features. Having input files that were almost, but not quite, like a programming language made me realize, however, that what I really wanted was a programming language--no matter how many features I added, I always wanted one more "simple" thing. As far as parallelization, I quickly realized that I had a problem: I needed a parallel FFT, and the only ones that were available were proprietary (non-portable), used incompatible data distribution formats, and were often designed to be called only from special languages like HPF. That, plus my dissatisfaction with the available free FFTs, led me to embark on a side project (with my friend Matteo Frigo of MIT's Laboratory for Computer Science) to develop a new, free FFT library, FFTW, that included portable parallel routines. Also, I decided that attempting to support too many models of parallel programming (threads, MPI, Cray shmem) in one program resulted in a mess; it was better to stick with MPI (supporting running only a single process too, of course). Another mistake I discovered was that I allowed the eigensolver to get too intertwined with the specific problem of Maxwell's equations--the eigensolver knew about the data structures for the fields, etcetera, making it difficult to plug in replacement eigensolvers, test things in isolation, or to implement features like the "targeted" eigensolver of Photonic-Bands (which diagonalizes a different operator). The whole program was too mired in complexity. Finally, in the interim I had learned about block eigensolver algorithms. Not only can such algorithms leverage prepackaged, highly-optimized routines like BLAS and LAPACK, but they also promised to be inherently more suited for parallelization (since they remove the serial process of solving for the bands one by one). All of these things convinced me that I needed to rewrite the code again from scratch.
So, I started work on the new package (in Spring 1998), this time determined to develop and debug each component (matrix operations, eigensolver, maxwell operator, user input) in isolation. At the same time, I was thinking about how I would implement the user control language, wanting to develop a general tool that could be applied to other problems and software in our research. It seemed clear that, in order to get other people to use it in their programs, as well as to avoid a lot of the manual labor that went into my previous effort, I wanted to automatically generate the user/program interface from an abstract specification. As for the control language, I briefly considered implementing my own, but was happily led to GNU Guile instead, which gave me a powerful language with little effort. So armed, I set out to write libctl (which generates the Guile interface from an abstract Scheme specification), partially as an experiment to see how hard it would be and what the result would look like. After a weekend of work, it was obvious that I had a powerful tool; I spent couple of weeks adding some finishing touches, writing documentation, and so on, and proudly showed it off to my groupmates in the hope that they could use it for their programs. Without a real example of a program using libctl, however, it was hard to convince them to plug a scripting language into their existing, working codes. So, I went back to puttering at my eigensolvers.
Of course, all this time I was allegedly doing real research, and long periods would go by with little progress on the Photonic-Bands package. The original, Fortran code was still working, and in time one learned to bear its quirks and limitations with stoicism, although we cringed every time we had to show it to anyone else. By the summer of 1999, I had a working block eigensolver (supporting several iteration-scheme variants), a Maxwell operator to plug into the eigensolver (including a "targeted" operator, whose convergence I was unhappy with), and a test program to do convergence experiments on bands of a Bragg mirror. I hadn't attached any general user interface, field output, or other necessary components. At this point, Dr. Doug Allan of Corning (a former student of Prof. Joannopoulos), heard about the new code--in particular, the targeted eigensolver--and began clamoring to try it out. Not put off by my excuses, he asked for a copy of my current code, regardless of its status, to play with. Not wanting to refuse, but aghast at the prospect of someone seeing my masterpiece only half-painted, I told him to give me a week...in which time I added the interface and discovered that I had a useful tool. Over the next week, I added many features, fixed bugs, and wrote documentation, drawing near to a release at last...
Installation
| |
| MIT Photonic Bands |
| Download |
| Release notes |
| MPB manual |
| Introduction |
| Installation |
| User Tutorial |
| Data Analysis Tutorial |
| User Reference |
| Developer Information |
| Acknowledgements |
| License and Copyright |
In this section, we outline the procedure for installing the MIT Photonic-Bands package. Mainly, this consists of downloading and installing various prerequisites. As much as possible, we have attempted to take advantage of existing packages such as BLAS, LAPACK, FFTW, and GNU Guile, in order to make our code smaller, more robust, faster, and more flexible. Unfortunately, this may make the installation of MPB more complicated if you do not already have these packages.
You will also need an ANSI C compiler, of course (gcc is fine), and installation will be easiest on a UNIX-like system (Linux is fine). In the following list, some of the packages are dependent upon packages listed earlier, so you should install them in more-or-less the order given.
Note: Many of these libraries may be available in precompiled binary form, especially for GNU/Linux systems. Be aware, however, that library binary packages often come in two parts, library and library-dev, and both are required to compile programs using it.
Note: It is important that you use the same Fortran compiler to compile Fortran libraries (like LAPACK) and for configuring MPB. Different Fortran compilers often have incompatible linking schemes. (The Fortran compiler for MPB can be set via the F77 environment variable.)
Unix Installation Basics
Installation Paths
First, let's review some important information about installing software on Unix systems, especially in regards to installing software in non-standard locations. (None of these issues are specific to this software, but they've caused a lot of confusion among our beloved users.)
Most of the software below, including this software, installs under /usr/local by default; that is, libraries go in /usr/local/lib, programs in /usr/local/bin, etc. If you don't have root privileges on your machine, you may need to install somewhere else, e.g. under $HOME/install (the install/ subdirectory of your home directory). Most of the programs below use a GNU-style configure script, which means that all you would do to install there would be:
./configure --prefix=$HOME/install
when configuring the program; the directories $HOME/install/lib etc. are created automatically as needed.
Paths for Configuring
There are two further complications. First, if you install in a non-standard location (and /usr/local is considered non-standard by some proprietary compilers), you will need to tell the compilers where to find the libraries and header files that you dutifully installed. You do this by setting two environment variables:
setenv LDFLAGS "-L/usr/local/lib" setenv CPPFLAGS "-I/usr/local/include"
Of course, substitute whatever installation directory you used. Do this before you run the configure scripts, etcetera. You may need to include multiple -L and -I flags (separated by spaces) if your machine has stuff installed in several non-standard locations. Bourne shell users (e.g. bash or ksh) should use the "export FOO=bar" syntax instead of csh's "setenv FOO bar", of course.
You might also need to update your PATH so that you can run the executables you installed (although /usr/local/bin is in the default PATH on many systems). e.g. if we installed in our home directory as described above, we would do:
setenv PATH "$HOME/install/bin:$PATH"
Paths for Running (Shared Libraries)
Second, many of the packages installed below (e.g. Guile) are installed as shared libraries. You need to make sure that your runtime linker knows where to find these shared libraries. The bad news is that every operating system does this in a slightly different way. The good news is that, when you run make install for the packages involving shared libraries, the output includes the necessary instructions specific to your system, so pay close attention! It will say something like "add LIBDIR to the <code>foobar environment variable</code>," where LIBDIR will be your library installation directory (e.g. /usr/local/lib) and foobar is some environment variable specific to your system (e.g. LD_LIBRARY_PATH on some systems, including Linux). For example, you might do:
setenv foobar "/usr/local/lib:$foobar"
Note that we just add to the library path variable, and don't replace it in case it contains stuff already. If you use Linux and have root privileges, you can instead simply run /sbin/ldconfig, first making sure that a line "/usr/local/lib" (or whatever) is in /etc/ld.so.conf.
If you don't want to type these commands every time you log in, you can put them in your ~/.cshrc file (or ~/.profile, or ~/.bash_profile, depending on your shell).
Fun with Fortran
Our software, along with many of the libraries it calls, is written in C or C++, but it also calls libraries such as BLAS and LAPACK (see below) that are usually compiled from Fortran. This can cause some added difficulty because of the various linking schemes used by Fortran compilers. Our configure script attempts to detect the Fortran linking scheme automatically, but in order for this to work you must use the same Fortran compiler and options with our software as were used to compile BLAS/LAPACK.
By default, our software looks for a vendor Fortran compiler first (f77, xlf, etcetera) and then looks for GNU g77. In order to manually specify a Fortran compiler foobar you would configure it with './configure F77=foobar ...'.
If, when you compiled BLAS/LAPACK, you used compiler options that alter the linking scheme (e.g. g77's -fcase-upper or -fno-underscoring), you will need to pass the same flags to our software via './configure FFLAGS="...flags..." ...'.
Picking a compiler
It is often important to be consistent about which compiler you employ. This is especially true for C++ software. To specify a particular C compiler foo, configure with ./configure CC=foo; to specify a particular C++ compiler foo++, configure with ./configure CXX=foo++; to specify a particular Fortran compiler foo90, configure with ./configure F77=foo90.
GNU/Linux and BSD binary packages
If you are installing on your personal GNU/Linux or BSD machine, then precompiled binary packages are likely to be available for many of these packages, and may even have been included with your system. On Debian systems, the packages are in .deb format and the built-in apt-get program can fetch them from a central repository. On Red Hat, SuSE, and most other Linux-based systems, binary packages are in RPM format and can be fetched (if they did not come with your system) from places like rpmfind.net. OpenBSD has its "ports" system, and so on.
Do not compile something from source if an official binary package is available. For one thing, you're just creating pain for yourself. Worse, the binary package may already be installed, in which case installing a different version from source will just cause trouble.
One thing to watch out for is that libraries like LAPACK, Guile, HDF5, etcetera, will often come split up into two (or more) packages: e.g. a guile package and a guile-devel package. You need to install both of these to compile software using the library.
BLAS and LAPACK
MPB requires the BLAS and LAPACK libraries for matrix computations.
BLAS
The first thing you must have on your system is a BLAS implementation. "BLAS" stands for "Basic Linear Algebra Subroutines," and is a standard interface for operations like matrix multiplication. It is designed as a building-block for other linear-algebra applications, and is used both directly by our code and in LAPACK (see below). By using it, we can take advantage of many highly-optimized implementations of these operations that have been written to the BLAS interface. (Note that you will need implementations of BLAS levels 1-3.)
You can find more BLAS information, as well as a basic implementation, on the BLAS Homepage. Once you get things working with the basic BLAS implementation, it might be a good idea to try and find a more optimized BLAS code for your hardware. Vendor-optimized BLAS implementations are available as part of the Intel MKL, HP CXML, IBM ESSL, SGI sgimath, and other libraries. Recently, there has also been work on self-optimizing BLAS implementations that can achieve performance competitive with vendor-tuned codes; see the ATLAS homepage. ATLAS works well, but it does take some time to compile.
Note that the generic BLAS does not come with a Makefile; compile it with something like:
mkdir blas && cd blas # the BLAS archive does not create its own directory get http://www.netlib.org/blas/blas.tgz gunzip blas.tgz tar xf blas.tar f77 -c -O3 *.f # compile all of the .f files to produce .o files ar rv libblas.a *.o # combine the .o files into a library su -c "cp libblas.a /usr/local/lib" # switch to root and install
(Replace -O3 with your favorite optimization options. On Linux, I use g77 -O3 -fomit-frame-pointer -funroll-loops, with -malign-double -mcpu=i686 on a Pentium II.) Note that MPB looks for the standard BLAS library with -lblas, so the library file should be called libblas.a and reside in a standard directory like /usr/local/lib. (See also below for the --with-blas=lib option to MPB's configure script, to manually specify a library location.)
LAPACK
LAPACK, the Linear Algebra PACKage, is a standard collection of routines, built on BLAS, for more-complicated (dense) linear algebra operations like matrix inversion and diagonalization. You can download LAPACK from the LAPACK Home Page.
Note that our software looks for LAPACK by linking with -llapack. This means that the library must be called liblapack.a and be installed in a standard directory like /usr/local/lib (alternatively, you can specify another directory via the LDFLAGS environment variable as described earlier). (See also below for the --with-lapack=lib option to our configure script, to manually specify a library location.)
MPI (for parallel MPB)
Optionally, MPB is able to run on a distributed-memory parallel machine, and to do this we use the standard MPI message-passing interface. You can learn about MPI from the MPI Home Page. Most commercial supercomputers already have an MPI implementation installed; two free MPI implementations that you can install yourself are MPICH and LAM; a promising new implementation is Open MPI. MPI is not required to compile the ordinary, uniprocessor version of our software.
In order for the MPI version of our software to run successfully, we have a slightly nonstandard requirement: each process must be able to read from the disk. (This way, Guile can boot for each process and they can all read your control file in parallel.) Many (most?) commercial supercomputers, Linux Beowulf clusters, etcetera, satisfy this requirement.
Also, in order to get good performance, you'll need fast interconnect hardware such as Myrinet; 100Mbps Ethernet probably won't cut it. The speed bottleneck comes from the FFT, which requires most of the communications in MPB. FFTW's MPI transforms (see below) come with benchmark programs that will give you a good idea of whether you can get speedups on your system. Of course, even with slow communications, you can still benefit from the memory savings per CPU for large problems.
As described below, when you configure MPB with MPI support (--with-mpi), it installs itself as mpb-mpi. See also the [user-ref.html#mpb-mpi user reference] section for information on using MPB on parallel machines.
MPI support in MPB is thanks to generous support from Clarendon Photonics.
HDF5 (optional, strongly advised)
We require a portable, standard binary format for outputting the electromagnetic fields and similar volumetric data, and for this we use HDF. (If you don't have HDF5, you can still compile MPB, but you won't be able to output the fields or the dielectric function.)
HDF is a widely-used, free, portable library and file format for multi-dimensional scientific data, developed in the National Center for Supercomputing Applications (NCSA) at the University of Illinois (UIUC, home of the Fighting Illini and alma mater to many of this author's fine relatives). You can get HDF and learn about it on the HDF Home Page.
There are two incompatible versions of HDF, HDF4 and HDF5 (no, not HDF1 and HDF2). We require the newer version, HDF5, which is supported by a number scientific of visualization tools, including our own h5utils utilities.
HDF5 includes parallel I/O support under MPI, which can be enabled by configuring it with --enable-parallel. (You may also have to set the CC environment variable to mpicc.) Unfortunately, the parallel HDF5 library then does not work with serial code, so you have may have to choose one or the other.
We have some hacks in our MPB and Meep software so that they can do parallel I/O even with the serial HDF5 library; these hacks work okay when you are using a small number of processors, but on large supercomputers we strongly recommend using the parallel HDF5.
Note: version 1.8 of HDF5 changed the programming interface, breaking older programs like MPB (versions 1.4.x and earlier) and Meep (versions 0.11.x and earlier). As a workaround, if you have HDF5 1.8 or later and would like to compile these versions of MPB and/or Meep, include -DH5_USE_16_API=1 in the CPPFLAGS (see below).
FFTW
FFTW is a self-optimizing, portable, high-performance FFT implementation, including both serial and parallel FFTs. You can download FFTW and find out more about it from the FFTW Home Page.
Important: the current release of MPB requires FFTW 2.x, and will not compile with FFTW 3.x.
If you want to use MPB on a parallel machine with MPI, you will also need to install the MPI FFTW libraries (this just means including --enable-mpi in the FFTW configure flags).
GNU Readline (optional)
GNU Readline is a library to provide command-line history, tab-completion, emacs keybindings, and other shell-like niceties to command-line programs. This is an optional package, but one that can be used by Guile (see below) if it is installed; we recommend getting it. You can download Readline from the GNU ftp site. (Readline is typically preinstalled on GNU/Linux systems).
GNU Guile
GNU Guile is an extension/scripting language implementation based on Scheme, and we use it to provide a rich, fully-programmable user interface with minimal effort. It's free, of course, and you can download it from the Guile Home Page. (Guile is typically included with GNU/Linux systems.)
- Important: Most Linux distributions come with Guile already installed (check by seeing whether you can run
guile --versionfrom the command line). In that case, do not install your own version of Guile from source — having two versions of Guile on the same system will cause problems. However, by default most distributions install only the Guile libraries and not the programming headers — to compile libctl and MPB, you should install the "guile-devel" or "guile-dev" package (e.g. with RPM) that came with your system.
GNU Autoconf (optional)
If you want to be a developer of the MPB package (as opposed to merely a user), you will also need the GNU Autoconf program. Autoconf is a portability tool that generates configure scripts to automatically detect the capabilities of a system and configure a package accordingly. You can find out more at the Autoconf Home Page (autoconf is typically installed by default on Linux systems). In order to install Autoconf, you will also need the GNU m4 program if you do not already have it (see the GNU m4 Home Page).
libctl
Instead of using Guile directly, we separated much of the user interface code into a package called libctl, in the hope that this might be more generally useful. libctl automatically handles the communication between the program and Guile, converting complicated data structures and so on, to make it even easier to use Guile to control scientific applications. Download libctl from the libctl page, unpack it, and run the usual configure, make, make install sequence. You'll also want to browse the libctl manual, as this will give you a general overview of what the user interface will be like.
If you are not the system administrator of your machine, and/or want to install libctl somewhere else (like your home directory), you can do so with the standard --prefix=dir option to configure (the default prefix is /usr/local). In this case, however, you'll need to specify the location of the libctl shared files for the MPB or Meep package, using the --with-libctl=dir/share/libctl option to our configure script.
MIT Photonic Bands
Okay, if you've made it all the way here, you're ready to install the MPB package and start cranking out eigenmodes. (You can download the latest version and read this manual at the MIT Photonic-Bands Homepage.) Once you've unpacked it, just run:
./configure make
to configure and compile the package (see below to install). Hopefully, the configure script will correctly detect the BLAS, FFTW, etcetera libraries that you've dutifully installed, as well as the C compiler and so on, and the make compilation will proceed without a hitch. If not, it's a Simple Matter of Programming to correct the problem. configure accepts several flags to help control its behavior. Some of these are standard, like --prefix=dir to specify and installation directory prefix, and some of them are specific to the MPB package (./configure --help for more info). The configure flags specific to MPB are:
-
--with-inv-symmetry - Assume [user-ref.html#inv-symmetry inversion symmetry] in the dielectric function, allowing us to use real fields (in Fourier space) instead of complex fields. This gives a factor of 2 benefit in speed and memory. In this case, the MPB program will be installed as
mpbiinstead ofmpb, so that you can have versions both with and without inversion symmetry installed at the same time. To install bothmpbandmpbi, you should do:
./configure make su -c "make install" make distclean ./configure --with-inv-symmetry make su -c "make install"
-
--with-hermitian-eps - Support the use of [user-ref.html#dielectric-anisotropic complex-hermitian dielectric tensors] (corresponding to magnetic materials, which break inversion symmetry).
-
--enable-single - Use single precision (C
float) instead of the default double precision (Cdouble) for computations. (Not recommended.) -
--without-hdf5 - Don't use the HDF5 library for field and dielectric function output. (In which case, no field output is possible.)
-
--with-mpi - Attempt to compile a [user-ref.html#mpb-mpi parallel version of MPB] using MPI; the resulting program will be installed as
mpb-mpi. Requires [#mpi MPI] and [#fftw MPI FFTW] libraries to be installed, as described above.
Does not compile the serial MPB, or mpb-data; if you want those, you should make distclean and compile/install them separately.
--with-mpi can be used along with --with-inv-symmetry, in which case the program is installed as mpbi-mpi (try typing that five times quickly).
-
--with-libctl=dir - If libctl was installed in a nonstandard location (i.e. neither
/usrnor/usr/local), you need to specify the location of the libctl directory,dir. This is eitherprefix/share/libctl, whereprefixis the installation prefix of libctl, or the original libctl source code directory. -
--with-blas=lib - The
configurescript automatically attempts to detect accelerated BLAS libraries, like DXML (DEC/Alpha), SCSL and SGIMATH (SGI/MIPS), ESSL (IBM/PowerPC), ATLAS, and PHiPACK. You can, however, force a specific library name to try via--with-blas=lib. -
--with-lapack=lib - Cause the
configurescript to look for a LAPACK library calledlib(the default is to use-llapack). -
--disable-checks - Disable runtime checks. (Not recommended; the disabled checks shouldn't take up a significant amount of time anyway.)
-
--enable-prof - Compile for performance profiling.
-
--enable-debug - Compile for debugging, adding extra runtime checks and so on.
-
--enable-debug-malloc - Use special memory-allocation routines for extra debugging (to check for array overwrites, memory leaks, etcetera).
-
--with-efence - More debugging: use the Electric Fence library, if available, for extra runtime array bounds-checking.
You can further control configure by setting various environment variables, such as:
-
CC: the C compiler command -
CFLAGS: the C compiler flags (defaults to-O3). -
CPPFLAGS:-Idirflags to tell the C compiler additional places to look for header files. -
LDFLAGS:-Ldirflags to tell the linker additional places to look for libraries. -
LIBS: additional libraries to link against.
Once compiled, the main program (as opposed to various test programs) resides in the mpb-ctl/ subdirectory, and is called mpb. You can install this program under /usr/local (or elsewhere, if you used the --prefix flag for configure), by running:
su -c "make install"
The "su" command is to switch to root for installation into system directories. You can just do make install if you are installing into your home directory instead.
If you make a mistake (e.g. you forget to specify a needed -Ldir flag) or in general want to start over from a clean slate, you can restore MPB to a pristine state by running:
make distclean
User Tutorial
| |
| MIT Photonic Bands |
| Download |
| Release notes |
| MPB manual |
| Introduction |
| Installation |
| User Tutorial |
| Data Analysis Tutorial |
| User Reference |
| Developer Information |
| Acknowledgements |
| License and Copyright |
In this section, we'll walk through the process of computing the band structure and outputting some fields for a two-dimensional photonic crystal using the MIT Photonic-Bands package. This should give you the basic idea of how it works and some of the things that are possible. Here, we tell you the truth, but not the whole truth. The MPB User Reference section gives a more complete, but less colloquial, description of the features supported by Photonic-Bands. See also the [analysis-tutorial.html data analysis tutorial] in the next section for more examples, focused on analyzing and visualizing the results of MPB calculations.
The ctl File
The use of the Photonic-Bands package revolves around the control file, abbreviated "ctl" and typically called something like foo.ctl (although you can use any filename extension you wish). The ctl file specifies the geometry you wish to study, the number of eigenvectors to compute, what to output, and everything else specific to your calculation. Rather than a flat, inflexible file format, however, the ctl file is actually written in a scripting language. This means that it can be everything from a simple sequence of commands setting the geometry, etcetera, to a full-fledged program with user input, loops, and anything else that you might need.
Don't worry, though—simple things are simple (you don't need to be a Real Programmer), and even there you will appreciate the flexibility that a scripting language gives you. (e.g. you can input things in any order, without regard for whitespace, insert comments where you please, omit things when reasonable defaults are available...)
The ctl file is actually implemented on top of the libctl library, a set of utilities that are in turn built on top of the Scheme language. Thus, there are three sources of possible commands and syntax for a ctl file:
- Scheme, a powerful and beautiful programming language developed at MIT, which has a particularly simple syntax: all statements are of the form
(function arguments...). We run Scheme under the GNU Guile interpreter (designed to be plugged into programs as a scripting and extension language). You don't need to learn much Scheme for a basic ctl file, but it is always there if you need it; you can learn more about it from these Guile and Scheme links. - libctl, a library that we built on top of Guile to simplify communication between Scheme and scientific computation software. libctl sets the basic tone of the user interface and defines a number of useful functions. See the libctl pages.
- The MIT Photonic-Bands package, which defines all the interface features that are specific to photonic band structure calculations. This manual is primarily focused on documenting these features.
It would be an excellent idea at this point for you to go read the libctl manual, particularly the [[[libctl Basic User Experience]] section, which will give you an overview of what the user interface is like, provide a crash course in the Scheme features that are most useful here, and describe some useful general features. We're not going to repeat this material (much), so learn it now!
Okay, let's continue with our tutorial. The MIT Photonic-Bands (MPB) program is normally invoked by running something like:
unix% mpb foo.ctl >& foo.out
which reads the ctl file foo.ctl and executes it, saving the output to the file foo.out. (Some sample ctl files are provided for you in the mpb-ctl/examples/ directory.) However, as you already know (since you obediently read the libctl manual, right?), if you invoke mpb with no arguments, you are dropped into an interactive mode in which you can type commands and see their results immediately. Why don't you do that right now, in your terminal window? Then, you can paste in the commands from the tutorial as you follow it and see what they do. Isn't this fun?
Our First Band Structure
As our beginning example, we'll compute the band structure of a two-dimensional square lattice of dielectric rods in air (see our online textbook, ch. 5). In our control file, we'll first specify the parameters and geometry of the simulation, and then tell it to run and give us the output.
All of the parameters (each of which corresponds to a Scheme variable) have default setting, so we only need to specify the ones we need to change. (For a complete listing of the parameter variables and their current values, along with some other info, type (help) at the guile> prompt.) One of the parameters, num-bands, controls how many bands (eigenstates) are computed at each k point. If you type num-bands at the prompt, it will return the current value, 1--this is too small; let's set it to a larger value:
(set! num-bands 8)
This is how we change the value of variables in Scheme (if you type num-bands now, it will return 8). The next thing we want to set (although the order really doesn't matter), is the set of k points (Bloch wavevectors) we want to compute the bands at. This is controlled by the variable k-points, a list of 3-vectors (which is initially empty). We'll set it to the corners of the irreducible Brillouin zone of a square lattice, Gamma, X, M, and Gamma again:
(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
Notice how we construct a list, and how we make 3-vectors; notice also how we can break things into multiple lines if we want, and that a semicolon (';') marks the start of a comment. Typically, we'll want to also compute the bands at a lot of intermediate k points, so that we see the continuous band structure. Instead of manually specifying all of these intermediate points, however, we can just call one of the functions provided by libctl to interpolate them for us:
(set! k-points (interpolate 4 k-points))
This takes the k-points and linearly interpolates four new points between each pair of consecutive points; if we type k-points now at the prompt, it will show us all 16 points in the new list:
(#(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))
As is described below, all spatial vectors in the program are specified in the basis of the lattice directions normalized to basis-size lengths (default is unit-normalized), while the k points are specified in the basis of the (unnormalized) reciprocal lattice vectors. In this case, we don't have to specify the lattice directions, because we are happy with the defaults--the lattice vectors default to the cartesian unit axes (i.e. the most common case, a square/cubic lattice). (The reciprocal lattice vectors in this case are also the unit axes.) We'll see how to change the lattice vectors in later subsections.
Now, we want to set the geometry of the system--we need to specify which objects are in the primitive cell of the lattice, centered on the origin. This is controlled by the variable geometry, which is a list of geometric objects. As you know from reading the libctl documentation, objects (more complicated, structured data types), are created by statements of the form (make type (property1 value1) (property2 value2) ...). There are various kinds (sub-classes) of geometric object: cylinders, spheres, blocks, ellipsoids, and perhaps others in the future. Right now, we want a square lattice of rods, so we put a single dielectric cylinder at the origin:
(set! geometry (list (make cylinder
(center 0 0 0) (radius 0.2) (height infinity)
(material (make dielectric (epsilon 12))))))
Here, we've set several properties of the cylinder: the center is the origin, its radius is 0.2, and its height (the length along its axis) is infinity. Another property, the material, it itself an object--we made it a dielectric with the property that its epsilon is 12. There is another property of the cylinder that we can set, the direction of its axis, but we're happy with the default value of pointing in the z direction.
All of the geometric objects are ostensibly three-dimensional, but since we're doing a two-dimensional simulation the only thing that matters is their intersection with the xy plane (z=0). Speaking of which, let us set the dimensionality of the system. Normally, we do this when we define the size of the computational cell, controlled by the geometry-lattice variable, an object of the lattice class: we can set some of the dimensions to have a size no-size, which reduces the dimensionality of the system.
(set! geometry-lattice (make lattice (size 1 1 no-size)))
Here, we define a 1x1 two-dimensional cell (defaulting to square). This cell is discretized according to the resolution variable, which defaults to 10 (pixels/lattice-unit). That's on the small side, and this is only a 2d calculation, so let's increase the resolution:
(set! resolution 32)
This results in a 32x32 computational grid. For efficient calculation, it is best to make the grid sizes a power of two, or factorizable into powers of small primes (such as 2, 3, 5 and 7). As a rule of thumb, you should use a resolution of at least 8 in order to obtain reasonable accuracy.
Now, we're done setting the parameters--there are other parameters, but we're happy with their default values for now. At this point, we're ready to go ahead and compute the band structure. The simplest way to do this is to type (run). Since this is a two-dimensional calculation, however, we would like to split the bands up into TE- and TM-polarized modes, and we do this by invoking (run-te) and (run-tm).
These produce a lot of output, showing you exactly what the code is doing as it runs. Most of this is self-explanatory, but we should point out one output in particular. Among the output, you should see lines like:
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
These lines are designed to allow you to easily extract the band-structure information and import it into a spreadsheet for graphing. They comprise the k point index, the k components and magnitude, and the frequencies of the bands, in comma-delimited format. Each line is prefixed by "tefreqs" for TE bands, "tmfreqs" for TM bands, and "freqs" for ordinary bands (produced by (run)), and using this prefix you can extract the data you want from the output by passing it through a program like grep. For example, if you had redirected the output to a file foo.out (as described earlier), you could extract the TM bands by running grep tmfreqs foo.out at the Unix prompt. Note that the output includes a header line, like:
tefreqs:, k index, kx, ky, kz, kmag/2pi, band 1, band 2, band 3, band 4, band 5, band 6, band 7, band 8
explaining each column of data. Another output of the run is the list of band gaps detected in the computed bands. For example the (run-tm) output includes the following gap output:
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%
This data is also stored in the variable gap-list, which is a list of (gap-percent gap-min gap-max) lists. It is important to realize, however, that this band-gap data may include "false positives," from two possible sources:
- If two bands cross, a false gap may result because the code computes the gap (connecting the dots in the band diagram) by assuming that bands never cross. Such false gaps are typically quite small (< 1%). To be sure of what's going on, you should either look at the symmetry of the modes involved or compute k points very close to the crossing (although even if the crossing occurs precisely at one of your k-points, there usually won't be an exact degeneracy for numerical reasons).
- One typically computes band diagrams by considering k-points around the boundary of the irreducible Brillouin zone. It is possible (though rare) that the band edges may occur at points in the interior of the Brillouin zone. To be absolutely sure you have a band gap (and of its size), you should compute the frequencies for points inside the Brillouin zone, too.
You've computed the band structure, and extracted the eigenfrequencies for each k point. But what if you want to see what the fields look like, or check that the dielectric function is what you expect? To do this, you need to output HDF files for these functions (HDF is a binary format for multi-dimensional scientific data, and can be read by many visualization programs). (We output files in HDF5 format, where the filenames are suffixed by ".h5".)
When you invoke one of the run functions, the dielectric function in the unit cell is automatically written to the file "epsilon.h5". To output the fields or other information, you need to pass one or more arguments to the run function. For example:
(run-tm output-efield-z) (run-te (output-at-kpoint (vector3 0.5 0 0) output-hfield-z output-dpwr))
This will output the electric (E) field z components for the TM bands at all k-points; and the magnetic (H) field z components and electric field energy densities (D power) for the TE bands at the X point only. The output filenames will be things like "e.k12.b03.z.te.h5", which stands for the z component (.z) of the TE (.te) electric field (e) for the third band (.b03) of the twelfth k point (.k12). Each HDF5 file can contain multiple datasets; in this case, it will contain the real and imaginary parts of the field (in datasets "z.r" and "z.i"), and in general it may include other data too (e.g. output-efield outputs all the components in one file). See also the [analysis-tutorial.html Data Analysis] tutorial.
There are several other output functions you can pass, described in the [user-ref.html user reference section], like output-dfield, output-hpwr, and output-dpwr-in-objects. Actually, though, you can pass in arbitrary functions that can do much more than just output the fields--you can perform arbitrary analyses of the bands (using functions that we will describe later).
Instead of calling one of the run functions, it is also possible to call lower-level functions of the code directly, to have a finer control of the computation; such functions are described in the reference section.
A Few Words on Units
In principle, you can use any units you want with the Photonic-Bands package. You can input all of your distances and coordinates in furlongs, for example, as long as you set the size of the lattice basis accordingly (the basis-size property of geometry-lattice). We do not recommend this practice, however.
Maxwell's equations possess an important property--they are scale-invariant (see our online textbook, ch. 2). If you multiply all of your sizes by 10, the solution scales are simply multiplied by 10 likewise (while the frequencies are divided by 10). So, you can solve a problem once and apply that solution to all length-scales. For this reason, we usually pick some fundamental lengthscale a of a structure, such as its lattice constant (unit of periodicity), and write all distances in terms of that. That is, we choose units so that a is unity. Then, to apply to any physical system, one simply scales all distances by a. This is what we have done in the preceding and following examples, and this is the default behavior of Photonic-Bands: the lattice constant is one, and all coordinates are scaled accordingly.
As has been mentioned already, (nearly) all 3-vectors in the program are specified in the basis of the lattice vectors normalized to lengths given by basis-size, defaulting to the unit-normalized lattice vectors. That is, each component is multiplied by the corresponding basis vector and summed to compute the corresponding cartesian vector. It is worth noting that a basis is not meaningful for scalar distances (such as the cylinder radius), however: these are just the ordinary cartesian distances in your chosen units of a.
Note also that the k-points, as mentioned above, are an exception: they are in the basis of the reciprocal lattice vectors (see our online textbook, appendix B). If a given dimension has size no-size, its reciprocal lattice vector is taken to be 2*pi/a.
We provide conversion functions to transform vectors between the various bases.
The frequency eigenvalues returned by the program are in units of c/a, where c is the speed of light and a is the unit of distance. (Thus, the corresponding vacuum wavelength is a over the frequency eigenvalue.)
Bands of a Triangular Lattice
As a second example, we'll compute the TM band structure of a triangular lattice of dielectric rods in air. To do this, we only need to change the lattice, controlled by the variable geometry-lattice. We'll set it so that the first two basis vectors (the properties basis1 and basis2) point 30 degrees above and below the x axis, instead of their default value of the x and y axes:
(set! geometry-lattice (make lattice (size 1 1 no-size)
(basis1 (/ (sqrt 3) 2) 0.5)
(basis2 (/ (sqrt 3) 2) -0.5)))
We don't specify basis3, keeping its default value of the z axis. Notice that Scheme supplies us all the ordinary arithmetic operators and functions, but they use prefix (Polish) notation, in Scheme fashion. The basis properties only specify the directions of the lattice basis vectors, and not their lengths--the lengths default to unity, which is fine here.
The irreducible Brillouin zone of a triangular lattice is different from that of a square lattice, so we'll need to modify the k-points list accordingly:
(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))
Note that these vectors are in the basis of the new reciprocal lattice vectors, which are different from before. (Notice also the Scheme shorthand (/ 3), which is the same as (/ 1 3) or 1/3.)
All of the other parameters (geometry, num-bands, and grid-size) can remain the same as in the previous subsection, so we can now call (run-tm) to compute the bands. As it turns out, this structure has an even larger TM gap than the square lattice:
Gap from band 1 (0.275065617068082) to band 2 (0.446289918847647), 47.4729292989213%
Maximizing the First TM Gap
Just for fun, we will now show you a more sophisticated example, utilizing the programming capabilities of Scheme: we will write a script to choose the cylinder radius that maximizes the first TM gap of the triangular lattice of rods from above. All of the Scheme syntax here won't be explained, but this should give you a flavor of what is possible.
First, we will write the function that want to maximize, a function that takes a dielectric constant and returns the size of the first TM gap. This function will change the geometry to reflect the new radius, run the calculation, and return the size of the first gap:
(define (first-tm-gap r)
(set! geometry (list (make cylinder
(center 0 0 0) (radius r) (height infinity)
(material (make dielectric (epsilon 12))))))
(run-tm)
(retrieve-gap 1)) ; return the gap from TM band 1 to TM band 2
We'll leave most of the other parameters the same as in the previous example, but we'll also change num-bands to 2, since we only need to compute the first two bands:
(set! num-bands 2)
In order to distinguish small differences in radius during the optimization, it might seem that we have to increase the grid resolution, slowing down the computation. Instead, we can simply increase the mesh resolution. This is the size of the mesh over which the dielectric constant is averaged at each grid point, and increasing the mesh size means that the average index better reflects small changes in the structure.
(set! mesh-size 7) ; increase from default value of 3
Now, we're ready to maximize our function first-tm-gap. We could write a loop to do this ourselves, but libctl provides a built-in function (maximize function tolerance arg-min arg-max) to do it for us (using Brent's algorithm). So, we just tell it to find the maximum, searching in the range of radii from 0.1 to 0.5, with a tolerance of 0.1:
(define result (maximize first-tm-gap 0.1 0.1 0.5)) (print "radius at maximum: " (max-arg result) "\n") (print "gap size at maximum: " (max-val result) "\n")
(print is a function defined by libctl to apply the built-in display function to zero or more arguments.) After five iterations, the output is:
radius at maximum: 0.176393202250021 gap size at maximum: 48.6252611051049
The tolerance of 0.1 that we specified means that the true maximum is within 0.1 * 0.176393202250021, or about 0.02, of the radius found here. It doesn't make much sense here to specify a lower tolerance, since the discretization of the grid means that the code can't accurately distinguish small differences in radius.
Before we continue, let's reset mesh-size to its default value:
(set! mesh-size 3) ; reset to default value of 3
A Complete 2D Gap with an Anisotropic Dielectric
As another example, one which does not require so much Scheme knowledge, let's construct a structure with a complete 2D gap (i.e., in both TE and TM polarizations), in a somewhat unusual way: using an dielectric structure. An anisotropic dielectric presents a different dielectric constant depending upon the direction of the electric field, and can be used in this case to make the TE and TM polarizations "see" different structures.
We already know that the triangular lattice of rods has a gap for TM light, but not for TE light. The dual structure, a triangular lattice of holes, has a gap for TE light but not for TM light (at least for the small radii we will consider). Using an anisotropic dielectric, we can make both of these structures simultaneously, with each polarization seeing the structure that gives it a gap.
As before, our geometry will consist of a single cylinder, this time with a radius of 0.3, but now it will have an epsilon of 12 (dielectric rod) for TM light and 1 (air hole) for TE light:
(set! geometry (list (make cylinder
(center 0 0 0) (radius 0.3) (height infinity)
(material (make dielectric-anisotropic
(epsilon-diag 1 1 12))))))
Here, epsilon-diag specifies the diagonal elements of the dielectric tensor. The off-diagonal elements (specified by epsilon-offdiag) default to zero and are only needed when the principal axes of the dielectric tensor are different from the cartesian xyz axes.
The background defaults to air, but we need to make it a dielectric (epsilon of 12) for the TE light, so that the cylinder forms a hole. This is controlled via the default-material variable:
(set! default-material (make dielectric-anisotropic (epsilon-diag 12 12 1)))
Finally, we'll increase the number of bands back to eight and run the computation:
(set! num-bands 8) (run) ; just use run, instead of run-te or run-tm, to find the complete gap
The result, as expected, is a complete band gap:
Gap from band 2 (0.223977612336924) to band 3 (0.274704473679751), 20.3443687933601%
(If we had computed the TM and TE bands separately, we would have found that the lower edge of the complete gap in this case comes from the TM gap, and the upper edge comes from the TE gap.)
Finding a Point-defect State
Here, we consider the problem of finding a point-defect state in our square lattice of rods. This is a state that is localized in a small region by creating a point defect in the crystal — e.g., by removing a single rod. The resulting mode will have a frequency within, and be confined by, the gap (see our online textbook, ch. 5).
To compute this, we need a supercell of bulk crystal, within which to put the defect — we will use a 5x5 cell of rods. To do this, we must first increase the size of the lattice by five, and then add all of the rods. We create the lattice by:
(set! geometry-lattice (make lattice (size 5 5 no-size)))
Here, we have used the default orthogonal basis, but have changed the size of the cell. To populate the cell, we could specify all 25 rods manually, but that would be tedious — and any time you find yourself doing something tedious in a programming language, you're making a mistake. A better approach would be to write a loop, but in fact this has already been done for you. Photonic-Bands provides a function, geometric-objects-lattice-duplicates, that duplicates a list of objects over the lattice:
(set! geometry (list (make cylinder
(center 0 0 0) (radius 0.2) (height infinity)
(material (make dielectric (epsilon 12))))))
(set! geometry (geometric-objects-lattice-duplicates geometry))
There, now the geometry list contains 25 rods — the originally geometry list (which contained one rod, remember) duplicated over the 5x5 lattice.
To remove a rod, we'll just add another rod in the center of the cell with a dielectric constant of 1. When objects overlap, the later object in the list takes precedence, so we have to put the new rod at the end of geometry:
(set! geometry (append geometry
(list (make cylinder (center 0 0 0)
(radius 0.2) (height infinity)
(material air)))))
Here, we've used the Scheme append function to combine two lists, and have also snuck in the predefined material type air (which has an epsilon of 1).
We'll be frugal and only use only 16 points per lattice unit (resulting in an 80x80 grid) instead of the 32 from before:
(set! resolution 16)
Only a single k point is needed for a point-defect calculation (which, for an infinite supercell, would be independent of k):
(set! k-points (list (vector3 0.5 0.5 0)))
Unfortunately, for a supercell the original bands are folded many times over (in this case, 25 times), so we need to compute many more bands to reach the same frequencies:
(set! num-bands 50)
At this point, we can call (run-tm) to solve for the TM bands. It will take several minutes to compute, so be patient. Recall that the gap for this structure was for the frequency range 0.2812 to 0.4174. The bands of the solution include exactly one state in this frequency range: band 25, with a frequency of 0.378166. This is exactly what we should expect--the lowest band was folded 25 times into the supercell Brillouin zone, but one of these states was pushed up into the gap by the defect.
We haven't yet output any of the fields, but we don't have to repeat the run to do so. The fields from the last k-point computation remain in memory and can continue to be accessed and analyzed. For example, to output the electric field z component of band 25, we just do:
(output-efield-z 25)
That's right, the output functions that we passed to (run) in the first example are just functions of the band index that are called on each band. We can do other computations too, like compute the fraction of the electric field energy near the defect cylinder (within a radius 1.0 of the origin):
(get-dfield 25) ; compute the D field for band 25
(compute-field-energy) ; compute the energy density from D
(print
"energy in cylinder: "
(compute-energy-in-objects (make cylinder (center 0 0 0)
(radius 1.0) (height infinity)
(material air)))
"\n")
The result is 0.624794702341156, or over 62% of the field energy in this localized region; the field decays exponentially into the bulk crystal. The full range of available functions is described in the [user-ref.html reference section], but the typical sequence is to first load a field with a get- function and then to call other functions to perform computations and transformations on it.
Note that the compute-energy-in-objects returns the energy fraction, but does not itself print this value. This is fine when you are running interactively, in which case Guile always displays the result of the last expression, but when running as part of a script you need to explicitly print the result as we have done above with the print function. The "\n" string is newline character (like in C), to put subsequent output on a separate line.
Instead of computing all those bands, we can instead take advantage of a special feature of the MIT Photonic-Bands software that allows you to compute the bands closest to a "target" frequency, rather than the bands with the lowest frequencies. One uses this feature by setting the target-freq variable to something other than zero (e.g. the mid-gap frequency). In order to get accurate results, it's currently also recommended that you decrease the tolerance variable, which controls when convergence is judged to have occurred, from its default value of 1e-7:
(set! num-bands 1) ; only need to compute a single band, now! (set! target-freq (/ (+ 0.2812 0.4174) 2)) (set! tolerance 1e-8)
Now, we just call (run-tm) as before. Convergence requires more iterations this time, both because we've decreased the tolerance and because of the nature of the eigenproblem that is now being solved, but only by about 3-4 times in this case. Since we now have to compute only a single band, however, we arrive at an answer much more quickly than before. The result, of course, is again the defect band, with a frequency of 0.378166.
Tuning the Point-defect Mode
As another example utilizing the programming power of Scheme, we will write a script to "tune" the defect mode to a particular frequency. Instead of forming a defect by simply removing a rod, we can decrease the radius or the dielectric constant of the defect rod, thereby changing the corresponding mode frequency. In this case, we'll vary the dielectric constant, and try to find a mode with a frequency of, say, 0.314159 (a random number).
We could write a loop to search for this epsilon, but instead we'll use a root-finding function provided by libctl, (find-root function tolerance arg-min arg-max), that will solve the problem for us using a quadratically-convergent algorithm (Ridder's method). First, we need to define a function that takes an epsilon for the center rod and returns the mode frequency minus 0.314159; this is the function we'll be finding the root of:
(define old-geometry geometry) ; save the 5x5 grid with a missing rod
(define (rootfun eps)
; add the cylinder of epsilon = eps to the old geometry:
(set! geometry (append old-geometry
(list (make cylinder (center 0 0 0)
(radius 0.2) (height infinity)
(material (make dielectric
(epsilon eps)))))))
(run-tm) ; solve for the mode (using the targeted solver)
(print "epsilon = " eps " gives freq. = " (list-ref freqs 0) "\n")
(- (list-ref freqs 0) 0.314159)) ; return 1st band freq. - 0.314159
Now, we can solve for epsilon (searching in the range 1 to 12, with a fractional tolerance of 0.01) by:
(define rooteps (find-root rootfun 0.01 1 12)) (print "root (value of epsilon) is at: " rooteps "\n")
The sequence of dielectric constants that it tries, along with the corresponding mode frequencies, is:
| epsilon | frequency |
|---|---|
| 1 | 0.378165893321125 |
| 12 | 0.283987088221692 |
| 6.5 | 0.302998920718043 |
| 5.14623274327171 | 0.317371748739314 |
| 5.82311637163586 | 0.309702408341706 |
| 5.41898003340128 | 0.314169110036439 |
| 5.62104820251857 | 0.311893530112625 |
The final answer that it returns is an epsilon of 5.41986120170136 Interestingly enough, the algorithm doesn't actually evaluate the function at the final point; you have to do so yourself if you want to find out how close it is to the root. (Ridder's method successively reduces the interval bracketing the root by alternating bisection and interpolation steps. At the end, it does one last interpolation to give you its best guess for the root location within the current interval.) If we go ahead and evaluate the band frequency at this dielectric constant (calling (rootfun rooteps)), we find that it is 0.314159008193209, matching our desired frequency to nearly eight decimal places after seven function evaluations! (Of course, the computation isn't really this accurate anyway, due to the finite discretization.)
A slight improvement can be made to the calculation above. Ordinarily, each time you call the (run-tm) function, the fields are initialized to random values. It would speed convergence somewhat to use the fields of the previous calculation as the starting point for the next calculation. We can do this by instead calling a lower-level function, (run-parity TM false). (The first parameter is the polarization to solve for, and the second tells it not to reset the fields if possible.)
emacs and ctl
It is useful to have emacs use its scheme-mode for editing ctl files, so that hitting tab indents nicely, and so on. emacs does this automatically for files ending with ".scm"; to do it for files ending with ".ctl" as well, add the following lines to your ~/.emacs file:
(if (assoc "\\.ctl" auto-mode-alist)
nil
(nconc auto-mode-alist '(("\\.ctl" . scheme-mode))))
(Incidentally, emacs scripts are written in "elisp," a language closely related to Scheme.)
Data Analysis Tutorial
| |
| MIT Photonic Bands |
| Download |
| Release notes |
| MPB manual |
| Introduction |
| Installation |
| User Tutorial |
| Data Analysis Tutorial |
| User Reference |
| Developer Information |
| Acknowledgements |
| License and Copyright |
In the previous section, we focused on how to perform a calculation in MPB. Now, we'll give a brief tutorial on what you might do with the results of the calculations, and in particular how you might visualize the results. We'll focus on two systems, one two-dimensional and one three-dimensional.
Triangular Lattice of Rods
First, we'll return to the two-dimensional triangular lattice of rods in air from the tutorial (see also our online textbook, ch. 5). The control file for this calculation, which can also be found in mpb-ctl/examples/tri-rods.ctl, will consist of:
The tri-rods.ctl control file
(set! num-bands 8)
(set! geometry-lattice (make lattice (size 1 1 no-size)
(basis1 (/ (sqrt 3) 2) 0.5)
(basis2 (/ (sqrt 3) 2) -0.5)))
(set! geometry (list (make cylinder
(center 0 0 0) (radius 0.2) (height infinity)
(material (make dielectric (epsilon 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
(set! k-points (interpolate 4 k-points))
(set! resolution 32)
(run-tm (output-at-kpoint (vector3 (/ -3) (/ 3) 0)
fix-efield-phase output-efield-z))
(run-te)
Notice that we're computing both TM and TE bands (where we expect a gap in the TM bands), and are outputting the z component of the electric field for the TM bands at the K point. (The fix-efield-phase will be explained below.)
Now, run the calculation, directing the output to a file, by entering the following command at the Unix prompt:
unix% mpb tri-rods.ctl >& tri-rods.out
It should finish after a minute or two.
The tri-rods dielectric function
In most cases, the first thing we'll want to do is to look at the dielectric function, to make sure that we specified the correct geometry. We can do this by looking at the epsilon.h5 output file.
The first thing that might come to mind would be to examine epsilon.h5 directly, say by converting it to a PNG image with h5topng (from my free h5utils package), magnifying it by 3:
unix% h5topng -S 3 epsilon.h5
epsilon.png) is shown at right, and it initially seems wrong! Why is the rod oval-shaped and not circular? Actually, the dielectric function is correct, but the image is distorted because the primitive cell of our lattice is a rhombus (with 60-degree acute angles). Since the output grid of MPB is defined over the non-orthogonal unit cell, while the image produced by h5topng (and most other plotting programs) is square, the image is skewed.
We can fix the image in a variety of ways, but the best way is probably to use the mpb-data utility included (and installed) with MPB. mpb-data allows us to rearrange the data into a rectangular cell (-r) with the same area/volume, expand the data to include multiple periods (-m periods), and change the resolution per unit distance in each direction to a fixed value (-n resolution). man mpb-data or run mpb-data -h for more options. In this case, we'll rectify the cell, expand it to three periods in each direction, and fix the resolution to 32 pixels per a:
unix% mpb-data -r -m 3 -n 32 epsilon.h5
It's important to use -n when you use -r, as otherwise the non-square unit cell output by -r will have a different density of grid points in each direction, and appear distorted. The output of mpb-data is by default an additional dataset within the input file, as we can see by running h5ls:
unix% h5ls epsilon.h5
data Dataset {32, 32}
data-new Dataset {96, 83}
description Dataset {SCALAR}
lattice\ copies Dataset {3}
lattice\ vectors Dataset {3, 3}
mpb-data is the one called data-new. We can examine it by running h5topng again, this time explicitly specifying the name of the dataset (and no longer magnifying):
unix% h5topng epsilon.h5:data-new
The new epsilon.png output image is shown at right. As you can see, the rods are now circular as desired, and they clearly form a triangular lattice.
Gaps and and band diagram for tri-rods
At this point, let's check for band gaps by picking out lines with the word "Gap" in them:
unix% grep Gap tri-rods.out Gap from band 1 (0.275065617068082) to band 2 (0.446289918847647), 47.4729292989213% Gap from band 3 (0.563582903703468) to band 4 (0.593059066215511), 5.0968516236891% Gap from band 4 (0.791161222813268) to band 5 (0.792042731370125), 0.111357548663006% Gap from band 5 (0.838730315053238) to band 6 (0.840305955160638), 0.187683867865441% Gap from band 6 (0.869285340346465) to band 7 (0.873496724070656), 0.483294361375001% Gap from band 4 (0.821658212109559) to band 5 (0.864454087942874), 5.07627823271133%
The first five gaps are for the TM bands (which we ran first), and the last gap is for the TE bands. Note, however that the < 1% gaps are probably false positives due to band crossings, as described in the user tutorial. There are no complete (overlapping TE/TM) gaps, and the largest gap is the 47% TM gap as expected (see our online textbook, appendix C). (To be absolutely sure of this and other band gaps, we would also check k-points within the interior of the Brillouin zone, but we'll omit that step here.)
Next, let's plot out the band structure. To do this, we'll first extract the TM and TE bands as comma-delimited text, which can then be imported and plotted in our favorite spreadsheet/plotting program.
unix% grep tmfreqs tri-rods.out > tri-rods.tm.dat unix% grep tefreqs tri-rods.out > tri-rods.te.dat
The TM and TE bands are both plotted below against the "k index" column of the data, with the special k-points labelled. TM bands are shown in blue (filled circles) with the gaps shaded light blue, while TE bands are shown in red (hollow circles) with the gaps shaded light red.
Note that we truncated the upper frequencies at a cutoff of 1.0 c/a. Although some of our bands go above that frequency, we didn't compute enough bands to fill in all of the states in that range. Besides, we only really care about the states around the gap(s), in most cases.
The source of the TM gap: examining the modes
Now, let's actually examine the electric-field distributions for some of the bands (which were saved at the K point, remember). Besides looking neat, the field patterns will tell us about the characters of the modes and provide some hints regarding the origin of the band gap.
As before, we'll run mpb-data on the field output files (named e.k11.b*.z.tm.h5), and then run h5topng to view the results:
unix% mpb-data -r -m 3 -n 32 e.k11.b*.z.tm.h5 unix% h5topng -C epsilon.h5:data-new -c bluered -Z -d z.r-new e.k11.b*.z.tm.h5
Here, we've used the -C option to superimpose (crude) black contours of the dielectric function over the fields, -c bluered to use a blue-white-red color table, -Z to center the color scale at zero (white), and -d to specify the dataset name for all of the files at once. man h5topng for more information. (There are plenty of data-visualization programs available if you want more sophisticated plotting capabilities than what h5topng offers, of course; you can use h5totxt to convert the data to a format suitable for import into e.g. spreadsheets.)
Note that the dataset name is z.r-new, which is the real part of the z component of the output of mpb-data. (Since these are TM fields, the z component is the only non-zero part of the electric field.) The real and imaginary parts of the fields correspond to what the fields look like at half-period intervals in time, and in general they are different. However, at K they are redundant, due to the inversion symmetry of that k-point (proof left as an exercise for the reader). Usually, looking at the real parts alone gives you a pretty good picture of the state, especially if you use fix-efield-phase (see below), which chooses the phase to maximize the field energy in the real part. Sometimes, though, you have to be careful: if the real part happens to be zero, what you'll see is essentially numerical noise and you should switch to the imaginary part.
The resulting field images are shown below:
| TM band 1 | TM band 2 | TM band 3 | TM band 4 | TM band 5 | TM band 6 | TM band 7 | TM band 8 |
|---|---|---|---|---|---|---|---|
|
|
|
|
|
|
|
|
|
Your images should look the same as the ones above. If we hadn't included fix-efield-phase before output-efield-z in the ctl file, on the other hand, yours would have differed slightly (e.g. by a sign or a lattice shift), because by default the phase is random.
When we look at the real parts of the fields, we are really looking at the fields of the modes at a particular instant in time (and the imaginary part is half a period later). The point in time (relative to the periodic oscillation of the state) is determined by the phase of the eigenstate. The fix-efield-phase band function picks a canonical phase for the eigenstate, giving us a deterministic picture.
We can see several things from these plots:
First, the origin of the band gap is apparent. The lowest band is concentrated within the dielectric rods in order to minimize its frequency. The next bands, in order to be orthogonal, are forced to have a node within the rods, imposing a large "kinetic energy" (and/or "potential energy") cost and hence a gap (see our online textbook, ch. 5). Successive bands have more and more complex nodal structures in order to maintain orthogonality. (The contrasting absence of a large TE gap has to do with boundary conditions. The perpendicular component of the displacement field must be continuous across the dielectric boundary, but the parallel component need not be.)
We can also see the deep impact of symmetry on the states. The K point has C3v symmetry (not quite the full C6v symmetry of the dielectric structure). This symmetry group has only one two-dimensional representation--that is what gives rise to the degenerate pairs of states (2/3, 4/5, and 7/8), all of which fall into this "p-like" category (where the states transform like two orthogonal dipole field patterns, essentially). The other two bands, 1 and 6, transform under the trivial "s-like" representation (with band 6 just a higher-order version of 1).
Diamond Lattice of Spheres
"Then were the entrances of this world made narrow, full of sorrow and travail: they are but few and evil, full of perils, and very painful." (Ezra 4:7)
Now, let us turn to a three-dimensional structure, a diamond lattice of dielectric spheres in air (see our online textbook, ch. 6). The basic techniques to compute and analyze the modes of this structure are the same as in two dimensions, but of course, everything becomes more complicated in 3d. It's harder to find a structure with a complete gap, the modes are no longer polarized, the computations are far bigger, and visualization is much more difficult, for starters.
(The band gap of the diamond structure was first identified in: K. M. Ho, C. T. Chan, and C. M. Soukoulis, "Existence of a photonic gap in periodic dielectric structures," Phys. Rev. Lett. 65, 3152 (1990).)
The control file for this calculation, which can also be found in mpb-ctl/examples/diamond.ctl, consists of:
Diamond control file
(set! geometry-lattice (make lattice
(basis-size (sqrt 0.5) (sqrt 0.5) (sqrt 0.5))