MPB User Reference
From AbInitio
Revision as of 03:09, 3 April 2014 (edit) Stevenj (Talk  contribs) (→Storing and combining multiple fields) ← Previous diff 
Current revision (03:12, 3 April 2014) (edit) Stevenj (Talk  contribs) (→Group velocities) 

Line 382:  Line 382:  
; <code>(computegroupvelocitycomponent ''direction'')</code>  ; <code>(computegroupvelocitycomponent ''direction'')</code>  
: Returns a list of the groupvelocity components (units of ''c'') in the given ''direction'', one for each band at the lastcomputed kpoint. ''direction'' is a vector in the reciprocallattice basis (like the kpoints); its length is ignored. (This has the advantage of being three times faster than <code>computegroupvelocities</code>.)  : Returns a list of the groupvelocity components (units of ''c'') in the given ''direction'', one for each band at the lastcomputed kpoint. ''direction'' is a vector in the reciprocallattice basis (like the kpoints); its length is ignored. (This has the advantage of being three times faster than <code>computegroupvelocities</code>.)  
+  ; <code>(compute1groupvelocity ''whichband'')</code>, <code>(compute1groupvelocitycomponent ''direction'' ''whichband'')</code>  
+  : As above, but returns the group velocity or component thereof only for band <code>whichband</code>.  
== Field manipulation ==  == Field manipulation == 
Current revision
MIT Photonic Bands 
Download 
Release notes 
MPB manual 
Introduction 
Installation 
User Tutorial 
Data Analysis Tutorial 
User Reference 
Developer Information 
Acknowledgements 
License and Copyright 
Here, we document the features exposed to the user by the MIT PhotonicBands package. We do not document the Scheme language or the functions provided by libctl (see also the libctl User Reference section of the libctl manual).
Contents 
Input Variables
These are global variables that you can set to control various parameters of the PhotonicBands computation. They are also listed (along with their current values) by the (help)
command. In brackets after each variable is the type of value that it should hold. (The classes, complex datatypes like geometricobject
, are described in a later subsection. The basic datatypes, like integer
, boolean
, cnumber
, and vector3
, are defined by libctl.)

geometry
[list ofgeometricobject
class]  Specifies the geometric objects making up the structure being simulated. When objects overlap, later objects in the list take precedence. Defaults to no objects (empty list).

defaultmaterial
[materialtype
class]  Holds the default material that is used for points not in any object of the geometry list. Defaults to air (epsilon of 1). See also
epsiloninputfile
, below. 
ensureperiodicity
[boolean
]  If true (the default), then geometric objects are treated as if they were shifted by all possible lattice vectors; i.e. they are made periodic in the lattice.

geometrylattice
[lattice
class]  Specifies the basis vectors and lattice size of the computational cell (which is centered on the origin of the coordinate system). These vectors form the basis for all other 3vectors in the geometry, and the lattice size determines the size of the primitive cell. If any dimension of the lattice size is the special value
nosize
, then the dimension of the lattice is reduced (i.e. it becomes two or onedimensional). (That is, the dielectric function becomes twodimensional; it is still, in principle, a three dimensional system, and the kpoint vectors can be threedimensional.) Generally, you should make anynosize
dimension(s) perpendicular to the others. Defaults to the orthogonal xyz vectors of unit length (i.e. a square/cubic lattice). 
resolution
[number
orvector3
]  Specifies the computational grid resolution, in pixels per lattice unit (a lattice unit is one basis vector in a given direction). If
resolution
is avector3
, then specifies a different resolution for each direction; otherwise the resolution is uniform. (The grid size is then the product of the lattice size and the resolution, rounded up to the next positive integer.) Defaults to10
. You can call(optimizegridsize!)
after setting theresolution
andgeometrylattice
to adjust the grid size for maximal performance (this rounds the grid size in each direction to the nearest integer with small factors, to improve FFT speed). 
gridsize
[vector3
]  Specifies the size of the discrete computational grid along each of the lattice directions. Deprecated: the preferred method is to use the
resolution
variable, above, in which case thegridsize
defaults tofalse
. To get the grid size you should instead use the(getgridsize)
function. 
dimensions
[integer
]  Explicitly specifies the dimensionality of the simulation; if the value is less than 3, the sizes of the extra dimensions in
gridsize
are ignored (assumed to be one). Defaults to 3. Deprecated: the preferred method is to setgeometrylattice
to have size nosize in any unwanted dimensions. 
kpoints
[list ofvector3
]  List of Bloch wavevectors to compute the bands at, expressed in the basis of the reciprocal lattice vectors. The reciprocal lattice vectors are defined as follows (see our online textbook, appendix B): Given the lattice vectors
R_{i}
(not the basis vectors), the reciprocal lattice vectorG_{j}
satisfiesR_{i} * G_{j} = 2*Pi*delta_{i,j}
, wheredelta_{i,j}
is the Kronecker delta (1 fori = j
and 0 otherwise). (R_{i}
for anynosize
dimensions is taken to be the corresponding basis vector.) Normally, the wavevectors should be in the first Brillouin zone (see below).kpoints
defaults to none (empty list). 
numbands
[integer
]  Number of bands (eigenvectors) to compute at each k point. Defaults to 1.

targetfreq
[number
]  If zero, the lowestfrequency
numbands
states are solved for at each k point (ordinary eigenproblem). If nonzero, solve for thenumbands
states whose frequencies have the smallest absolute difference withtargetfreq
(special, "targeted" eigenproblem). Beware that the targeted solver converges more slowly than the ordinary eigensolver and may require a lowertolerance
to get reliable results. Defaults to 0. 
tolerance
[number
]  Specifies when convergence of the eigensolver is judged to have been reached (when the eigenvalues have a fractional change less than
tolerance
between iterations). Defaults to 1.0e7. 
filenameprefix
[string
]  A string prepended to all output filenames. Defaults to
"FILE"
, where your control file is FILE.ctl. You can change this tofalse
to use no prefix. 
epsiloninputfile
[string
]  If this string is not
""
(the default), then it should be the name of an HDF5 file whose first/only dataset defines a dielectric function (over some discrete grid). This dielectric function is then used in place ofdefaultmaterial
(i.e. where there are nogeometry
objects). The grid of the epsilon file dataset need not matchgridsize
; it is scaled and/or linearly interpolated as needed. The lattice vectors for the epsilon file are assumed to be the same asgeometrylattice
. [ Note that, even if the grid sizes match and there are no geometric objects, the dielectric function used by MPB will not be exactly the dielectric function of the epsilon file, unless you also setmeshsize
to1
(see above). ] 
eigensolverblocksize
[integer
]  The eigensolver uses a "block" algorithm, which means that it solves for several bands simultaneously at each kpoint.
eigensolverblocksize
specifies this number of bands to solve for at a time; if it is zero or >=numbands
, then all the bands are solved for at once. Ifeigensolverblocksize
is a negative number, n, then MPB will try to use nearlyequal blocksizes close to n. Making the block size a small number can reduce the memory requirements of MPB, but block sizes > 1 are usually more efficient (there is typically some optimum size for any given problem). Defaults to 11 (i.e. solve for around 11 bands at a time). 
simplepreconditioner?
[boolean
]  Whether or not to use a simplified preconditioner; defaults to
false
(this is fastest most of the time). (Turning this on increases the number of iterations, but decreases the time for each iteration.) 
deterministic?
[boolean
]  Since the fields are initialized to random values at the start of each run, there are normally slight differences in the number of iterations, etcetera, between runs. Setting
deterministic?
totrue
makes things deterministic (the default isfalse
). 
eigensolverflags
[integer
]  This variable is undocumented and reserved for use by Jedi Masters only.
Predefined Variables
Variables predefined for your convenience and amusement.

air
,vacuum
[materialtype
class]  Two aliases for a predefined material type with a dielectric constant of 1.

nothing
[materialtype
class]  A material that, effectively, punches a hole through other objects to the background (
defaultmaterial
orepsiloninputfile
). 
infinity
[number
]  A big number (1.0e20) to use for "infinite" dimensions of objects.
Output Variables
Global variables whose values are set upon completion of the eigensolver.

freqs
[list ofnumber
]  A list of the frequencies of each band computed for the last k point. Guaranteed to be sorted in increasing order. The frequency of band
b
can be retrieved via(listref freqs ( b 1))
. 
iterations
[integer
]  The number of iterations required for convergence of the last k point.

parity
[string
]  A string describing the current required parity/polarization (
"te"
,"zeven"
, etcetera, or""
for none), useful for prefixing output lines for grepping.
Yet more global variables are set by the run
function (and its variants), for use after run
completes or by a band function (which is called for each band during the execution of run
.

currentk
[vector3
]  The k point (from the
kpoints
list) most recently solved. 
gaplist
[list of(percent freqmin freqmax)
lists]  This is a list of the gaps found by the eigensolver, and is set by the
run
functions when two or more kpoints are solved. (It is the empty list if no gaps are found.) 
bandrangedata
[list of((min . kpoint) . (max . kpoint))
pairs (of pairs)]  For each band, this list contains the minimum and maximum frequencies of the band, and the associated k points where the extrema are achieved. Note that the bands are defined by sorting the frequencies in increasing order, so this can be confused if two bands cross.
Classes
Classes are complex datatypes with various "properties" which may have default values. Classes can be "subclasses" of other classes; subclasses inherit all the properties of their superclass, and can be used any place the superclass is expected. An object of a class is constructed with:
(make class (prop1 val1) (prop2 val2) ...)
See also the libctl manual.
MIT PhotonicBands defines several types of classes, the most numerous of which are the various geometric object classes. You can also get a list of the available classes, along with their property types and default values, at runtime with the (help)
command.
lattice
The lattice class is normally used only for the geometrylattice
variable, and specifies the three lattice directions of the crystal and the lengths of the corresponding lattice vectors.

lattice
 Properties:

basis1
,basis2
,basis3
[vector3
]  The three lattice directions of the crystal, specified in the cartesian basis. The lengths of these vectors are ignoredonly their directions matter. The lengths are determined by the
basissize
property, below. These vectors are then used as a basis for all other 3vectors in the ctl file. They default to the x, y, and z directions, respectively. 
basissize
[vector3
]  The components of
basissize
are the lengths of the three basis vectors, respectively. They default to unit lengths. 
size
[vector3
]  The size of the lattice (i.e. the length of the lattice vectors
R_{i}
, in which the crystal is periodic) in units of the basis vectors. Thus, the actual lengths of the lattice vectors are given by the components ofsize
multiplied by the components ofbasissize
. (Alternatively, you can think ofsize
as the vector between opposite corners of the primitive cell, specified in the lattice basis.) Defaults to unit lengths.
If any dimension has the special size nosize
, then the dimensionality of the problem is reduced by one; strictly speaking, the dielectric function is taken to be uniform along that dimension. (In this case, the nosize
dimension should generally be orthogonal to the other dimensions.)
materialtype
This class is used to specify the materials that geometric objects are made of. Currently, there are three subclasses, dielectric
, dielectricanisotropic
, and materialfunction
.

dielectric
 A uniform, isotropic, linear dielectric material, with one property:

epsilon
[number
]  The dielectric constant (must be positive). No default value. You can also use
(index n)
as a synonym for(epsilon (* n n))
.


dielectricanisotropic
 A uniform, possibly anisotropic, linear dielectric material. For this material type, you specify the (realsymmetric, or possibly complexhermitian) dielectric tensor (relative to the cartesian xyz axes):
This allows your dielectric to have different dielectric constants for fields polarized in different directions. The epsilon tensor must be positivedefinite (have all positive eigenvalues); if it is not, MPB exits with an error. (This does not imply that all of the entries of the epsilon matrix need be positive.) The components of the tensor are specified via three properties:

epsilondiag
[vector3
]  The diagonal elements (a b c) of the dielectric tensor. No default value.

epsilonoffdiag
[cvector3
]  The offdiagonal elements (u v w) of the dielectric tensor. Defaults to zero. This is a
cvector3
, which simply means that the components may be complex numbers (e.g.3+0.1i
). If nonzero imaginary parts are specified, then the dielectric tensor is complexhermitian. This is only supported when MPB is configured with thewithhermitianeps
flag. This is not dissipative (the eigenvalues of epsilon are real), but rather breaks timereversal symmetry, corresponding to a gyrotropic (magnetooptic) material (see our online textbook, ch. 2). Note that inversion symmetry may not mean what you expect for complexhermitian epsilon, so be cautious about usingmpbi
in this case. 
epsilonoffdiagimag
[vector3
]  Deprecated: The imaginary parts of the offdiagonal elements (u v w) of the dielectric tensor; defaults to zero. Setting the imaginary parts directly by specifying complex numbers in
epsilonoffdiag
is preferred.

For example, a material with a dielectric constant of 3.0 for TE fields (polarized in the xy plane) and 5.0 for TM fields (polarized in the z direction) would be specified via (make (dielectricanisotropic (epsilondiag 3 3 5)))
. Please be aware that not all 2d anisotropic dielectric structures will have TE and TM modes, however.

materialfunction
 This material type allows you to specify the material as an arbitrary function of position. (For an example of this, see the
braggsine.ctl
file in theexamples/
directory.) It has one property:
materialfunc
[function
]  A function of one argument, the position
vector3
(in lattice coordinates), that returns the material at that point. Note that the function you supply can return any material; wild and crazy users could even return anothermaterialfunction
object (which would then have its function invoked in turn).  Instead of
materialfunc
, you can useepsilonfunc
: forepsilonfunc
, you give it a function of position that returns the dielectric constant at that point.

Normally, the dielectric constant is required to be positive (or positivedefinite, for a tensor). However, MPB does have a somewhat experimental feature allowing negative dielectrics (e.g. in a plasma). To use it, call the function (allownegativeepsilon)
before (run)
. In this case, it will output the (real) frequency squared in place of the (possibly imaginary) frequencies. (Convergence will be somewhat slower because the eigenoperator is not positive definite.)
geometricobject
This class, and its descendants, are used to specify the solid geometric objects that form the dielectric structure being simulated. The base class is:

geometricobject
 Properties:

material
[materialtype
class]  The material that the object is made of (usually some sort of dielectric). No default value (must be specified).

center
[vector3
]  Center point of the object. No default value.

One normally does not create objects of type geometricobject
directly, however; instead, you use one of the following subclasses. Recall that subclasses inherit the properties of their superclass, so these subclasses automatically have the material
and center
properties (which must be specified, since they have no default values).
Recall that all 3vectors, including the center of an object, its axes, and so on, are specified in the basis of the normalized lattice vectors normalized to basissize
. Note also that 3vector properties can be specified by either (property (vector3 x y z))
or, equivalently, (property x y z)
.
In a twodimensional calculation, only the intersections of the objects with the xy plane are considered.

sphere
 A sphere. Properties:

radius
[number
]  Radius of the sphere. No default value.


cylinder
 A cylinder, with circular crosssection and finite height. Properties:

radius
[number
]  Radius of the cylinder's crosssection. No default value.

height
[number
]  Length of the cylinder along its axis. No default value.

axis
[vector3
]  Direction of the cylinder's axis; the length of this vector is ignored. Defaults to point parallel to the z axis.


cone
 A cone, or possibly a truncated cone. This is actually a subclass of
cylinder
, and inherits all of the same properties, with one additional property. The radius of the base of the cone is given by theradius
property inherited fromcylinder
, while the radius of the tip is given by the new property,radius2
. (Thecenter
of a cone is halfway between the two circular ends.)
radius2
[number
]  Radius of the tip of the cone (i.e. the end of the cone pointed to by the
axis
vector). Defaults to zero (a "sharp" cone).


block
 A parallelepiped (i.e., a brick, possibly with nonorthogonal axes). Properties:

size
[vector3
]  The lengths of the block edges along each of its three axes. Not really a 3vector (at least, not in the lattice basis), but it has three components, each of which should be nonzero. No default value.

e1
,e2
,e3
[vector3
]  The directions of the axes of the block; the lengths of these vectors are ignored. Must be linearly independent. They default to the three lattice directions.


ellipsoid
 An ellipsoid. This is actually a subclass of
block
, and inherits all the same properties, but defines an ellipsoid inscribed inside the block.
Here are some examples of geometric objects created using the above classes, assuming that the lattice directions (the basis) are just the ordinary unit axes, and m
is some material we have defined:
; A cylinder of infinite radius and height 0.25 pointing along the x axis, ; centered at the origin: (make cylinder (center 0 0 0) (material m) (radius infinity) (height 0.25) (axis 1 0 0))
; An ellipsoid with its long axis pointing along (1,1,1), centered on ; the origin (the other two axes are orthogonal and have equal ; semiaxis lengths): (make ellipsoid (center 0 0 0) (material m) (size 0.8 0.2 0.2) (e1 1 1 1) (e2 0 1 1) (e3 2 1 1))
; A unit cube of material m with a spherical air hole of radius 0.2 at ; its center, the whole thing centered at (1,2,3): (set! geometry (list (make block (center 1 2 3) (material m) (size 1 1 1)) (make sphere (center 1 2 3) (material air) (radius 0.2))))
Functions
Here, we describe the functions that are defined by the PhotonicBands package. There are many types of functions defined, ranging from utility functions for duplicating geometric objects to run functions that start the computation.
See also the reference section of the libctl manual, which describes a number of useful functions defined by libctl.
Geometry utilities
Some utility functions are provided to help you manipulate geometric objects:

(shiftgeometricobject obj shiftvector)
 Translate
obj
by the 3vectorshiftvector
. 
(geometricobjectduplicates shiftvector minmultiple maxmultiple obj)
 Return a list of duplicates of
obj
, shifted by various multiples ofshiftvector
(fromminmultiple
tomaxmultiple
, inclusive, in steps of 1). 
(geometricobjectsduplicates shiftvector minmultiple maxmultiple objlist)
 Same as
geometricobjectduplicates
, except operates on a list of objects,objlist
. If A appears before B in the input list, then all the duplicates of A appear before all the duplicates of B in the output list. 
(geometricobjectslatticeduplicates objlist [ ux uy uz ])
 Duplicates the objects in
objlist
by multiples of the lattice basis vectors, making all possible shifts of the "primitive cell" (see below) that fit inside the lattice cell. (This is useful for supercell calculations; see the [usertutorial.html tutorial].) The primitive cell to duplicate isux
byuy
byuz
, in units of the basis vectors. These three parameters are optional; any that you do not specify are assumed to be1
. 
(pointinobject? point obj)
 Returns whether or not the given 3vector
point
is inside the geometric objectobj
. 
(pointinperiodicobject? point obj)
 As
pointinobject?
, but also checks translations of the given object by the lattice vectors. 
(displaygeometricobjectinfo indentby obj)
 Outputs some information about the given
obj
, indented byindentby
spaces.
Coordinate conversion functions
The following functions allow you to easily convert back and forth between the lattice, cartesian, and reciprocal bases. (See also the note on units in the tutorial.)

(lattice>cartesian x)
,(cartesian>lattice x)
 Convert
x
between the lattice basis (the basis of the lattice vectors normalized tobasissize
) and the ordinary cartesian basis, wherex
is either avector3
or amatrix3x3
, returning the transformed vector/matrix. In the case of a matrix argument, the matrix is treated as an operator on vectors in the given basis, and is transformed into the same operator on vectors in the new basis. 
(reciprocal>cartesian x)
,(cartesian>reciprocal x)
 Like the above, except that they convert to/from reciprocal space (the basis of the reciprocal lattice vectors). Also, the cartesian vectors output/input are in units of 2π.

(reciprocal>lattice x)
,(lattice>reciprocal x)
 Convert between the reciprocal and lattice bases, where the conversion again leaves out the factor of 2π (i.e. the latticebasis vectors are assumed to be in units of 2π).
Also, a couple of rotation functions are defined, for convenience, so that you don't have to explicitly convert to cartesian coordinates in order to use libctl's rotatevector3
function (see the Libctl User Reference):

(rotatelatticevector3 axis theta v)
,(rotatereciprocalvector3 axis theta v)
 Like
rotatevector3
, except thataxis
andv
are specified in the lattice/reciprocal bases.
Usually, kpoints are specified in the first Brillouin zone, but sometimes it is convenient to specify an arbitrary kpoint. However, the accuracy of MPB degrades as you move farther from the first Brillouin zone (due to the choice of a fixed planewave set for a basis). This is easily fixed: simply transform the kpoint to a corresponding point in the first Brillouin zone, and a completely equivalent solution (identical frequency, fields, etcetera) is obtained with maximum accuracy. The following function accomplishes this:

(firstbrillouinzone k)
 Given a kpoint
k
(in the basis of the reciprocal lattice vectors, as usual), return an equivalent point in the first Brillouin zone of the current lattice (geometrylattice
).
Note that firstbrillouinzone
can be applied to the entire kpoints
list with the Scheme expression: (map firstbrillouinzone kpoints)
.
Run functions
These are functions to help you run and control the simulation. The ones you will most commonly use are the run
function and its variants. The syntax of these functions, and one lowerlevel function, is:

(run bandfunc ...)
 This runs the simulation described by the input parameters (see above), with no constraints on the polarization of the solution. That is, it reads the input parameters, initializes the simulation, and solves for the requested eigenstates of each kpoint. The dielectric function is outputted to "
epsilon.h5
" before any eigenstates are computed.run
takes as arguments zero or more "band functions"bandfunc
. A band function should be a function of one integer argument, the band index, so that(bandfunc whichband)
performs some operation on the bandwhichband
(e.g. outputting fields). After every kpoint, each band function is called for the indices of all the bands that were computed. Alternatively, a band function may be a "thunk" (function of zero arguments), in which case(bandfunc)
is called exactly once per kpoint. 
(runzeven bandfunc ...), (runzodd bandfunc ...)
 These are the same as the
run
function except that they constrain their solutions to have even and odd symmetry with respect to the z=0 plane. You should use these functions only for structures that are symmetric through the z=0 mirror plane, where the third basis vector is in the z direction (0,0,1) and is orthogonal to the other two basis vectors, and when the k vectors are in the xy plane. Under these conditions, the eigenmodes always have either even or odd symmetry. In two dimensions, even/odd parities are equivalent to TE/TM polarizations, respectively (and are often strongly analogous even in 3d). Such a symmetry classification is useful for structures such as waveguides and photoniccrystal slabs. (For example, see the paper by S. G. Johnson et al., "Guided modes in photonic crystal slabs," PRB 60, 5751, August 1999.) 
(runte bandfunc ...), (runtm bandfunc ...)
 These are the same as the
run
function except that they constrain their solutions to be TE and TMpolarized, respectively, in two dimensions. The TE and TM polarizations are defined has having electric and magnetic fields in the xy plane, respectively. Equivalently, the H/E field of TE/TM light has only a z component (making it easier to visualize).
These functions are actually equivalent to calling runzeven
and runzodd
, respectively.
Note that for the modes to be segregated into TE and TM polarizations, the dielectric function must have mirror symmetry for reflections through the xy plane. If you use anisotropic dielectrics, you should be aware that they break this symmetry if the z direction is not one of the principle axes. If you use runte
or runtm
in such a case of broken symmetry, MPB will exit with an error.

(runyeven bandfunc ...), (runyodd bandfunc ...)
 These functions are analogous to
runzeven
andrunzodd
, except that they constrain their solutions to have even and odd symmetry with respect to the y=0 plane. You should use these functions only for structures that are symmetric through the y=0 mirror plane, where the second basis vector is in the y direction (0,1,0) and is orthogonal to the other two basis vectors, and when the k vectors are in the xz plane. 
runyevenzeven
,runyevenzodd
,runyoddzeven
,runyoddzodd
,runteyeven
,runteyodd
,runtmyeven
,runtmyodd
 These
run
like functions combine theyeven
/yodd
constraints withzeven
/zodd
orte
/tm
. See alsorunparity
, below. 
(runparity p resetfields bandfunc ...)
 Like the
run
function, except that it takes two extra parameters, a parityp
and a boolean (true
/false
) valueresetfields
.p
specifies a parity constraint, and should be one of the predefined variables:
NOPARITY
: equivalent torun

EVENZ
(orTE
): equivalent torunzeven
orrunte

ODDZ
(orTM
): equivalent torunzodd
orruntm

EVENY
(likeEVENZ
but for y=0 plane) 
ODDY
(likeODDZ
but for y=0 plane)

It is possible to specify more than one symmetry constraint simultaneously by adding them, e.g. (+ EVENZ ODDY)
requires the fields to be even through z=0 and odd through y=0. It is an error to specify incompatible constraints (e.g. (+ EVENZ ODDZ)
). Important: if you specify the z/y parity, the dielectric structure (and the k vector) must be symmetric about the z/y=0 plane, respectively.
If resetfields
is false
, the fields from any previous calculation will be reused as the starting point from this calculation, if possible; otherwise, the fields are reset to random values. The ordinary run
functions use a default resetfields
oftrue
. Alternatively, resetfields
may be a string, the name of an HDF5 file to load the initial fields from (as exported by saveeigenvectors
, below).

(displayeigensolverstats)
 Display some statistics on the eigensolver convergence; this function is useful mainly for MPB developers in tuning the eigensolver.
Several band functions for outputting the eigenfields are defined for your convenience, and are described in the Band output functions section, below. You can also define your own band functions, and for this purpose the functions described in the section Field manipulation functions, below, are useful. A band function takes the form:
(define (mybandfunc whichband) ...do stuff here with band index whichband... )
Note that the output variable freqs
may be used to retrieve the frequency of the band (see above). Also, a global variable currentk
is defined holding the current kpoint vector from the kpoints
list.
There are also some even lowerlevel functions that you can call, although you should not need to do most of the time:

(initparams p resetfields?)
 Read the input variables and initialize the simulation in preparation for computing the eigenvalues. The parameters are the same as the first two parameters of
runparity
. This function must be called before any of the other simulation functions below. (Note, however, that therun
functions all callinitparams
.) 
(setparity p)
 After calling
initparams
, you can change the parity constraint without resetting the other parameters by calling this function. Beware that this does not randomize the fields (see below); you don't want to try to solve for, say, the TM eigenstates when the fields are initialized to TE states from a previous calculation. 
(randomizefields)
 Initialize the fields to random values.

(solvekpoint k)
 Solve for the requested eigenstates at the Bloch wavevector
k
.
The inverse problem: k as a function of frequency
MPB's (run)
function(s) and its underlying algorithms compute the frequency w
as a function of wavevector k
. Sometimes, however, it is desirable to solve the inverse problem, for k
at a given frequency w
. This is useful, for example, when studying coupling in a waveguide between different bands at the same frequency (frequency is conserved even when wavevector is not). One also uses k(w)
to construct wavevector diagrams, which aid in understanding diffraction (e.g. negativediffraction materials and superprisms). To solve such problems, therefore, we provide the findk
function described below, which inverts w(k)
via a few iterations of Newton's method (using the group velocity dw/dk
). Because it employs a rootfinding method, you need to specify bounds on k
and a crude initial guess (order of magnitude is usually good enough).

(findk p omega bandmin bandmax kdir tol kmagguess kmagmin kmagmax [bandfunc...])
 Find the wavevectors in the current geometry/structure for the bands from
bandmin
tobandmax
at the frequencyomega
along thekdir
direction in kspace. Returns a list of the wavevector magnitudes for each band; the actual wavevectors are(vector3scale magnitude (unitvector3 kdir))
. The arguments offindk
are:
p
: parity (same as first argument torunparity
, above). 
omega
: the frequency at which to find the bands 
bandmin
,bandmax
: the range of bands to solve for the wavevectors of (inclusive). 
kdir
: the direction in kspace in which to find the wavevectors. (The magnitude ofkdir
is ignored.) 
tol
: the fractional tolerance with which to solve for the wavevector;1e4
is usually sufficient. (Like thetolerance
input variable, this is only the tolerance of the numerical iteration...it does not have anything to do with e.g. the error from finite gridresolution
.) 
kmagguess
: an initial guess for the k magnitude (alongkdir
) of the wavevector atomega
. Can either be a list (one guess for each band frombandmin
tobandmax
) or a single number (same guess for all bands, which is usually sufficient). 
kmagmin
,kmagmax
: a range of k magnitudes to search; should be large enough to include the correct k values for all bands. 
bandfunc
: zero or more band functions, just as in(run)
, which are evaluated at the computed k points for each band.

The findk
routine also prints a line suitable for grepping:
kvals: omega, bandmin, bandmax, kdir1, kdir2, kdir3, k magnitudes...
Band/output functions
All of these are functions that, given a band index, output the corresponding field or compute some function thereof (in the primitive cell of the lattice). They are designed to be passed as band functions to the run
routines, although they can also be called directly. See also the section on field normalizations.

(outputhfield whichband)

(outputhfieldx whichband)

(outputhfieldy whichband)

(outputhfieldz whichband)
 Output the magnetic (H) field for
whichband
; either all or one of the (Cartesian) components, respectively. 
(outputdfield whichband)

(outputdfieldx whichband)

(outputdfieldy whichband)

(outputdfieldz whichband)
 Output the electric displacement (D) field for
whichband
; either all or one of the (Cartesian) components, respectively. 
(outputefield whichband)

(outputefieldx whichband)

(outputefieldy whichband)

(outputefieldz whichband)
 Output the electric (E) field for
whichband
; either all or one of the (Cartesian) components, respectively. 
(outputhpwr whichband)
 Output the timeaveraged magneticfield energy density (hpwr = for
whichband
. 
(outputdpwr whichband)
 Output the timeaveraged electricfield energy density (dpwr = ) for
whichband
. 
(fixhfieldphase whichband)

(fixdfieldphase whichband)

(fixefieldphase whichband)
 Fix the phase of the given eigenstate in a canonical way based on the given spatial field (see also
fixfieldphase
, below). Otherwise, the phase is random; these functions also maximize the real part of the given field so that one can hopefully just visualize the real part. To fix the phase for output, pass one of these functions torun
before the corresponding output function, e.g.(runtm fixdfieldphase outputdfieldz)
Although we try to maximize the "realness" of the field, this has a couple of limitations. First, the phase of the different field components cannot, of course, be chosen independently, so an individual field component may still be imaginary. Second, if you use mpbi
to take advantage of inversion symmetry in your problem, the phase is mostly determined elsewhere in the program; fix_fieldphase
in that case only determines the sign.
See also below for the outputpoynting
and outputtotpwr
functions to output the Poynting vector and the total electromagnetic energy density, respectively, and the outputchargedensity
function to output the bound charge density.
Sometimes, you only want to output certain bands. For example, here is a function that, given an band/output function like the ones above, returns a new output function that only calls the first function for bands with a large fraction of their energy in an object(s). (This is useful for picking out defect states in supercell calculations.)

(outputdpwrinobjects bandfunc minenergy objects...)
 Given a band function
bandfunc
, returns a new band function that only callsbandfunc
for bands having a fraction of their electricfield energy greater thanminenergy
inside the given objects (zero or more geometric objects). Also, for each band, prints the fraction of their energy in the objects in the following form (suitable for grepping):
dpwr:, bandindex, frequency, energyinobjects
outputdpwrinobjects
only takes a single band function as a parameter, but if you want it to call several band functions, you can easily combine them into one with the following routine:

(combinebandfunctions bandfuncs...)
 Given zero or more band functions, returns a new band function that calls all of them in sequence. (When passed zero parameters, returns a band function that does nothing.)
It is also often useful to output the fields only at a certain kpoint, to let you look at typical field patterns for a given band while avoiding gratuitous numbers of output files. This can be accomplished via:

(outputatkpoint kpoint bandfuncs...)
 Given zero or more band functions, returns a new band function that calls all of them in sequence, but only at the specified
kpoint
. For other kpoints, does nothing.
Miscellaneous functions

(retrievegap lowerband)
 Return the frequency gap from the band #
lowerband
to the band #(lowerband
+1), as a percentage of midgap frequency. The "gap" may be negative if the maximum of the lower band is higher than the minimum of the upper band. (The gap is computed from thebandrangedata
of the previous run.)
Parity
Given a set of eigenstates at a kpoint, MPB can compute their parities with respect to the z=0 or y=0 plane. The z/y parity of a state is defined as the expectation value (under the usual inner product) of the mirrorflip operation through z/y=0, respectively. For true even and odd eigenstates (see e.g. runzeven
and runzodd
), this will be +1 and 1, respectively; for other states it will be something in between.
This is useful e.g. when you have a nearly symmetric structure, such as a waveguide with a substrate underneath, and you want to tell which bands are evenlike (parity > 0) and oddlike (parity < 0). Indeed, any state can be decomposed into purely even and odd functions, with absolutevaluesquared amplitudes of (1+parity)/2 and (1parity)/2, respectively.

displayzparities
,displayyparities
 These are band functions, designed to be passed to
(run)
, which output all of the z/y parities, respectively, at each kpoint (in commadelimited format suitable for grepping). 
(computezparities)
 Returns a list of the parities about the z=0 plane, one number for each band computed at the last kpoint.

(computeyparities)
 Returns a list of the parities about the y=0 plane, one number for each band computed at the last kpoint.
(The reader should recall that the magnetic field is only a pseudovector, and is therefore multiplied by 1 under mirrorflip operations. For this reason, the magnetic field appears to have opposite symmetry from the electric field, but is really the same.)
Group velocities
Given a set of eigenstates at a given kpoint, MPB can compute their group velocities (the derivative of frequency with respect to wavevector) using the HellmanFeynmann theorem. Three functions are provided for this purpose, and we document them here from highestlevel to lowestlevel.

displaygroupvelocities
 This is a band function, designed to be passed to
(run)
, which outputs all of the group velocity vectors (in the Cartesian basis, in units of c) at each kpoint. 
(computegroupvelocities)
 Returns a list of groupvelocity vectors (in the Cartesian basis, units of c) for the bands at the lastcomputed kpoint.

(computegroupvelocitycomponent direction)
 Returns a list of the groupvelocity components (units of c) in the given direction, one for each band at the lastcomputed kpoint. direction is a vector in the reciprocallattice basis (like the kpoints); its length is ignored. (This has the advantage of being three times faster than
computegroupvelocities
.) 
(compute1groupvelocity whichband)
,(compute1groupvelocitycomponent direction whichband)
 As above, but returns the group velocity or component thereof only for band
whichband
.
Field manipulation
The PhotonicBands package provides a number of ways to take the field of a band and manipulate, process, or output it. These methods usually work in two stages. First, one loads a field into memory (computing it in position space) by calling one of the get
functions below. Then, other functions can be called to transform or manipulate the field.
The simplest class of operations involve only the currentlyloaded field, which we describe in the second subsection below. To perform more sophisticated operations, involving more than one field, one must copy or transform the current field into a new field variable, and then call one of the functions that operate on multiple field variables (described in the third subsection).
Field normalization
In order to perform useful operations on the fields, it is important to understand how they are normalized. We normalize the fields in the way that is most convenient for perturbation and coupledmode theory [c.f. SGJ et al., PRE 65, 066611 (2002)] (see also our online textbook, ch. 2), so that their energy densities have unit integral. In particular, we normalize the electric (), displacement () and magnetic () fields, so that:
where the integrals are over the computational cell. Note the volume element (the volume of a grid pixel/voxel). If you simply sum over all the grid points, therefore, you will get (# grid points) / (volume of cell).
Note that we have dropped the pesky factors of 1/2, π, etcetera from the energy densities, since these do not appear in e.g. perturbation theory, and the fields have arbitrary units anyway. The functions to compute/output energy densities below similarly use and without any prefactors.
Loading and manipulating the current field
In order to load a field into memory, call one of the get
functions follow. They should only be called after the eigensolver has run (or after initparams
, in the case of getepsilon
). One normally calls them after run
, or in one of the band functions passed to run
.

(gethfield whichband)
 Loads the magnetic (H) field for the band
whichband
. 
(getdfield whichband)
 Loads the electric displacement (D) field for the band
whichband
. 
(getefield whichband)
 Loads the electric (E) field for the band
whichband
. (This function actually callsgetdfield
followed bygetefieldfromdfield
, below.) 
(getchargedensity whichband)
 Loads the bound charge density for the band
whichband
. 
(getepsilon)
 Loads the dielectric function.
Once loaded, the field can be transformed into another field or a scalar field:

(getefieldfromdfield)
 Multiplies by the inverse dielectric tensor to compute the electric field from the displacement field. Only works if a D field has been loaded.

(fixfieldphase)
 Fix the currentlyloaded eigenstate's phase (which is normally random) in a canonical way, based on the spatial field (H, D, or E) that has currently been loaded. The phase is fixed to make the real part of the spatial field as big as possible (so that you can hopefully visualize just the real part of the field), and a canonical sign is chosen. See also the
fix*fieldphase
band functions, above, which are convenient wrappers aroundfixfieldphase

(computefieldenergy)
 Given the H or D fields, computes the corresponding energy density function (normalized by the total energy in H or D, respectively). Also prints the fraction of the field in each of its cartesian components in the following form (suitable for grepping):
fenergycomponents:, kindex, bandindex, xfraction, yfraction, zfraction
where f
is either h
or d
. The return value of computefieldenergy
is a list of 7 numbers: (U xr xi yr yi zr zi)
. U
is the total, unnormalized energy, which is in arbitrary units deriving from the normalization of the eigenstate (e.g. the total energy for H is always 1.0). xr
is the fraction of the energy in the real part of the field's x component, xi
is the fraction in the imaginary part of the x component, etcetera (yr + yi = yfraction
, and so on).

(computefielddivergence)
 Given a vector field, compute its divergence.
Various integrals and other information about the eigenstate can be accessed by the following functions, useful e.g. for perturbation theory. Functions dealing with the field vectors require a field to be loaded, and functions dealing with the energy density require an energy density to be loaded via computefieldenergy
.

(computeenergyindielectric mineps maxeps)
 Returns the fraction of the energy that resides in dielectrics with epsilon in the range
mineps
tomaxeps
. 
(computeenergyinobjects objects...)
 Returns the fraction of the energy inside zero or more geometric objects.

(computeenergyintegral f)

f
is a function(f u eps r)
that returns a number given three parameters:u
, the energy density at a point;eps
, the dielectric constant at the same point; andr
, the position vector (in lattice coordinates) of the point.computeenergyintegral
returns the integral off
over the unit cell. (The integral is computed simply as the sum over the grid points times the volume of a grid pixel/voxel.) This can be useful e.g. for perturbationtheory calculations. 
(computefieldintegral f)
 Like
computeenergyintegral
, butf
is a function(f F eps r)
that returns a number (possibly complex) whereF
is the complex field vector at the given point. 
(getepsilonpoint r)
 Given a position vector
r
(in lattice coordinates), return the interpolated dielectric constant at that point. (Since MPB uses a an effective dielectric tensor internally, this actually returns the mean dielectric constant.) 
(getepsiloninversetensorpoint r)
 Given a position vector
r
(in lattice coordinates), return the interpolated inverse dielectric tensor (a 3x3 matrix) at that point. (Near a dielectric interface, the effective dielectric constant is a tensor even if you input only scalar dielectrics; see the epsilon overview for more information.) The returned matrix may be complexHermetian if you are employing magnetic materials. 
(getenergypoint r)
 Given a position vector
r
(in lattice coordinates), return the interpolated energy density at that point. 
(getfieldpoint r)
 Given a position vector
r
(in lattice coordinates), return the interpolated (complex) field vector at that point. 
(getblochfieldpoint r)
 Given a position vector
r
(in lattice coordinates), return the interpolated (complex) Bloch field vector at that point (this is the field without the exp(ikx) envelope).
Finally, we have the following functions to output fields (either the vector fields, the scalar energy density, or epsilon), with the option of outputting several periods of the lattice.

(outputfield [ nx [ ny [ nz ] ] ])

(outputfieldx [ nx [ ny [ nz ] ] ])

(outputfieldy [ nx [ ny [ nz ] ] ])

(outputfieldz [ nx [ ny [ nz ] ] ])
 Output the currentlyloaded field. The optional (as indicated by the brackets) parameters
nx
,ny
, andnz
indicate the number of periods to be outputted along each of the three lattice directions. Omitted parameters are assumed to be 1. For vector fields,outputfield
outputs all of the (Cartesian) components, while the other variants output only one component. 
(outputepsilon [ nx [ ny [ nz ] ] ])
 A shortcut for calling
getepsilon
followed byoutputfield
. Note that, because epsilon is a tensor, a number of datasets are outputted in"epsilon.h5"
:
"data"
: 3/trace(1/epsilon) 
"epsilon.{xx,xy,xz,yy,yz,zz}"
: the (Cartesian) components of the (symmetric) dielectric tensor. 
"epsilon_inverse.{xx,xy,xz,yy,yz,zz}"
: the (Cartesian) components of the (symmetric) inverse dielectric tensor.

Storing and combining multiple fields
In order to perform operations involving multiple fields, e.g. computing the Poynting vector , they must be stored in field variables. Field variables come in three flavors, realscalar (rscalar) fields, complexscalar (cscalar) fields, and complexvector (cvector) fields. There is a predefined field variable curfield
representing the currentlyloaded field (see above), and you can "clone" it to create more field variables with one of:

(fieldmake f)
 Return a new field variable of the same type and size as the field variable
f
. Does not copy the field contents (seefieldcopy
andfieldset!
, below). 
(rscalarfieldmake f)

(cscalarfieldmake f)

(cvectorfieldmake f)
 Like
fieldmake
, but return a realscalar, complexscalar, or complexvector field variable, respectively, of the same size asf
but ignoringf
's type. 
(cvectorfieldnonbloch! f)
 By default, complex vector fields are assumed to be Blochperiodic (and are multiplied by e^{ikx} in output routines). This function tells MPB that the complex vector field f should never be multiplied by Bloch phases.

(fieldset! fdest fsrc)
 Set
fdest
to store the same field values asfsrc
, which must be of the same size and type. 
(fieldcopy f)
 Return a new field variable that is exact copy of
f
; this is equivalent to callingfieldmake
followed byfieldset!
. 
(fieldload f)
 Loads the field
f
as the current field, at which point you can use all of the functions in the previous section to operate on it or output it.
Once you have stored the fields in variables, you probably want to compute something with them. This can be done in three ways: combining fields into new fields with fieldmap!
(e.g. combine E and H to ), integrating some function of the fields with integratefields
(e.g. to compute coupling integrals for perturbation theory), and getting the field values at arbitrary points with *fieldgetpoint
(e.g. to do a line or surface integral). These three functions are described below:

(fieldmap! fdest func [f1 f2 ...])
 Compute the new field
fdest
to be(func f1val f2val ...)
at each point in the grid, wheref1val
etcetera is the corresponding value off1
etcetera. All the fields must be of the same size, and the argument and return types offunc
must match those of thef1...
andfdest
fields, respectively.fdest
may be the same field as one of thef1...
arguments. Note: all fields are without Bloch phase factors exp(ikx). 
(integratefields func [f1 f2 ...])
 Compute the integral of the function
(func r [f1 f2 ...])
over the computational cell, wherer
is the position (in the usual lattice basis) andf1
etc. are fields (which must all be of the same size). (The integral is computed simply as the sum over the grid points times the volume of a grid pixel/voxel.) Note: all fields are without Bloch phase factors exp(ikx). See also the note below. 
(cvectorfieldgetpoint f r)

(cvectorfieldgetpointbloch f r)

(rscalarfieldgetpoint f r)
 Given a position vector
r
(in lattice coordinates), return the interpolated field cvector/rscalar fromf
at that point.cvectorfieldgetpointbloch
returns the field without the exp(ikx) Bloch wavevector, in analogue togetblochfieldpoint
.
You may be wondering how to get rid of the field variables once you are done with them: you don't, since they are garbage collected automatically.
We also provide functions, in analogue to e.g getefield
and outputefield
above, to "get" various useful functions as the current field and to output them to a file:

(getpoynting whichband)
 Loads the Poynting vector for the band
whichband
, the flux density of electromagnetic energy flow, as the current field. 1/2 of the real part of this vector is the timeaverage flux density (which can be combined with the imaginary part to determine the amplitude and phase of the timedependent flux). 
(outputpoynting whichband)

(outputpoyntingx whichband)

(outputpoyntingy whichband)

(outputpoyntingz whichband)
 Output the Poynting vector field for
whichband
; either all or one of the (Cartesian) components, respectively. 
(gettotpwr whichband)
 Load the timeaveraged electromagneticfield energy density () for
whichband
. (If you multiply the real part of the Poynting vector by a factor of 1/2, above, you should multiply by a factor of 1/4 here for consistency.) 
(outputtotpwr whichband)
 Output the timeaveraged electromagneticfield energy density (above) for
whichband
. 
(outputchargedensity whichband)
 Output the bound charge density (above) for
whichband
.
As an example, below is the Scheme source code for the getpoynting
function, illustrating the use of the various field functions:
(define (getpoynting whichband) (getefield whichband) ; put E in curfield (let ((e (fieldcopy curfield))) ; ... and copy to local var. (gethfield whichband) ; put H in curfield (fieldmap! curfield ; write ExH to curfield (lambda (e h) (vector3cross (vector3conj e) h)) e curfield) (cvectorfieldnonbloch! curfield))) ; see below
Stored fields and Bloch phases
Complex vector fields like E and H as computed by MPB are physically of the Bloch form: exp(ikx) times a periodic function. What MPB actually stores, however, is just the periodic function, the Bloch envelope, and only multiplies by exp(ikx) for when the fields are output or passed to the user (e.g. in integration functions). This is mostly transparent, with a few exceptions noted above for functions that do not include the exp(ikx) Bloch phase (it is somewhat faster to operate without including the phase).
On some occasions, however, when you create a field with fieldmap!
, the resulting field should not have any Bloch phase. For example, for the Poynting vector E^{*}xH, the exp(ikx) cancels because of the complex conjugation. After creating this sort of field, we must use the special function cvectorfieldnonbloch!
to tell MPB that the field is purely periodic:

(cvectorfieldnonbloch! f)
 Specify that the field
f
is not of the Bloch form, but rather that it is purely periodic.
Currently, all fields must be either Bloch or nonBloch (i.e. periodic), which covers most physically meaningful possibilities.
There is another wrinkle: even for fields in Bloch form, the exp(ikx) phase currently always uses the current kpoint, even if the field was computed from another kpoint. So, if you are performing computations combining fields from different kpoints, you should take care to always use the periodic envelope of the field, putting the Bloch phase in manually if necessary.
Manipulating the raw eigenvectors
MPB also includes a few lowlevel routines to manipulate the raw eigenvectors that it computes in a transverse planewave basis.
The most basic operations involve copying, saving, and restoring the current set of eigenvectors or some subset thereof:

(geteigenvectors firstband numbands)
 Return an eigenvector object that is a copy of
numbands
current eigenvectors starting atfirstband
. e.g. to get a copy of all of the eigenvectors, use(geteigenvectors 1 numbands)
. 
(seteigenvectors ev firstband)
 Set the current eigenvectors, starting at
firstband
, to those in theev
eigenvector object (as returned bygeteigenvectors
). (Does not work if the grid sizes don't match) 
(loadeigenvectors filename)

(saveeigenvectors filename)
 Read/write the current eigenvectors (raw planewave amplitudes) to/from an HDF5 file named
filename
. Instead of usingloadeigenvectors
directly, you can pass thefilename
as theresetfields
parameter ofrunparity
, above. (Loaded eigenvectors must be of the same size (same grid size and #bands) as the current settings.)
Currently, there's only one other interesting thing you can do with the raw eigenvectors, and that is to compute the dotproduct matrix between a set of saved eigenvectors and the current eigenvectors. This can be used, e.g., to detect band crossings or to set phases consistently at different k points. The dot product is returned as a "sqmatrix" object, whose elements can be read with the sqmatrixsize
and sqmatrixref
routines.

(doteigenvectors ev firstband)
 Returns a sqmatrix object containing the dot product of the saved eigenvectors
ev
with the current eigenvectors, starting atfirstband
. That is, the (i,j)th output matrix element contains the dot product of the (i+1)th vector ofev
(conjugated) with the (firstband
+j)th eigenvector. Note that the eigenvectors, when computed, are orthonormal, so the dot product of the eigenvectors with themselves is the identity matrix. 
(sqmatrixsize sm)
 Return the size n of an nxn sqmatrix
sm
. 
(sqmatrixref sm i j)
 Return the (
i
,j
)th element of the nxn sqmatrixsm
, where {i
,j
} range from 0..n1.
Inversion Symmetry
If you configure
MPB with the withinvsymmetry
flag, then the program is configured to assume inversion symmetry in the dielectric function. This allows it to run at least twice as fast and use half as much memory as the more general case. This version of MPB is by default installed as mpbi
, so that it can coexist with the usual mpb
program.
Inversion symmetry means that if you transform (x,y,z) to (x,y,z) in the coordinate system, the dielectric structure is not affected. Or, more technically, that (see our online textbook, ch. 3):
where the conjugation is significant for complexhermitian dielectric tensors. This symmetry is very common; all of the examples in this manual have inversion symmetry, for example.
Note that inversion symmetry is defined with respect to a specific origin, so that you may "break" the symmetry if you define a given structure in the wrong way—this will prevent mpbi
from working properly. For example, the diamond structure that we considered earlier would not have possessed inversion symmetry had we positioned one of the "atoms" to lie at the origin.
You might wonder what happens if you pass a structure lacking inversion symmetry to mpbi
. As it turns out, mpbi
only looks at half of the structure, and infers the other half by the inversion symmetry, so the resulting structure always has inversion symmetry, even if its original description did not. So, you should be careful, and look at the epsilon.h5
output to make sure it is what you expected.
Parallel MPB
We provide two methods by which you can parallelize MPB. The first, using MPI, is the most sophisticated and potentially provides the greatest and most general benefits. The second, which involves a simple script to split e.g. the kpoints
list among several processes, is less general but may be useful in many cases.
MPB with MPI parallelization
If you configure
MPB with the withmpi
flag, then the program is compiled to take advantage of distributedmemory parallel machines with MPI, and is installed as mpbmpi
. (See also the installation section.) This means that computations will (potentially) run more quickly and take up less memory per processor than for the serial code. (Normally, you should also install the serial version of MPB, if only to get the mpbdata
program, which is not installed with mpbmpi
.)
Using the parallel MPB is almost identical to using the serial version(s), with a couple of minor exceptions. The same ctl files should work for both. Running a program that uses MPI requires slightly different invocations on different systems, but will typically be something like:
unix% mpirun np 4 mpbmpi foo.ctl
to run on e.g. 4 processors. A second difference is that 1D systems are currently not supported in the MPI code, but the serial code should be fast enough for those anyway. A third difference is that the output HDF5 files (epsilon, fields, etcetera) from mpbmpi
have their first two dimensions (x and y) transposed; i.e. they are output as YxXxZ arrays. This doesn't prevent you from visualizing them, but the coordinate system is lefthanded; to untranspose the data, you can process it with mpbdata
and the T
option (in addition to any other options).
In order to get optimal benefit (time and memory savings) from mpbmpi
, the first two dimensions (n_{x} and n_{y}) of your grid should both be divisible by the number of processes. If you violate this constraint, MPB will still work, but the load balance between processors will be uneven. At worst, e.g. if either n_{x} or n_{y} is smaller than the number of processes, then some of the processors will be idle for part (or all) of the computation. When using inversion symmetry (mpbimpi
) for 2D grids only, the optimal case is somewhat more complicated: n_{x} and (n_{y}/2 + 1), not n_{y}, should both be divisible by the number of processes.
mpbmpi
divides each band at each kpoint between the available processors. This means that, even if you have only a single kpoint (e.g. in a defect calculation) and/or a single band, it can benefit from parallelization. Moreover, memory usage per processor is inversely proportional to the number of processors used. For sufficiently large problems, the speedup is also nearly linear.
MPI support in MPB is thanks to generous support from Clarendon Photonics.
Alternative parallelization: mpbsplit
There is an alternative method of parallelization when you have multiple k points: do each kpoint on a different processor. This does not provide any memory benefits, and does not allow one kpoint to benefit by starting with the fields of the previous kpoint, but is easy and may be the only effective way to parallelize calculations for small problems. This method also does not require MPI: it can utilize the unmodified serial mpb
program. To make it even easier, we supply a simple script called mpbsplit
(or mpbisplit
) to break the kpoints
list into chunks for you. Running:
unix% mpbsplit numsplit foo.ctl
will break the kpoints
list in foo.ctl
into numsplit
moreorless equal chunks, launch numsplit
processes of mpb
in parallel to process each chunk, and output the results of each in order. (Each process is an ordinary mpb
execution, except that it numbers its kpoints
depending upon which chunk it is in, so that output files will not overwrite one another and you can still grep
for frequencies as usual.)
Of course, this will only benefit you on a system where different processes will run on different processors, such as an SMP or a cluster with automatic process migration (e.g. MOSIX). mpbsplit
is actually a trivial shell script, though, so you can easily modify it if you need to use a special command to launch processes on other processors/machines (e.g. via GNU Queue).
The general syntax for mpbsplit
is:
unix% mpbsplit numsplit mpbarguments...
where all of the arguments following numsplit
are passed along to mpb
. What mpbsplit
technically does is to set the MPB variable ksplitnum
to numsplit
and ksplitindex
to the index (starting with 0) of the chunk for each process. If you want, you can use these variables to divide the problem in some other way and then reset them to 1
and , respectively.