Input file syntax

Input file

The structure to be simulated as well as all simulation parameters have to be specified in an XML file which can be edited by the users. The contents of such a file, e.g. THz_QCL_Fathololoumi_OpticsExpress2012.xml, are described in the following.

General syntax

Quotation marks are necessary if an equal sign (=) is present, e.g. AlloyX="Barrier".

<?xml version="1.0" encoding="utf-8"?>

<nextnano.MSB Version="1.0.1">  <!-- This number must be consistent to the version of the executable used. -->

Comment Section

<Header>                                       <!-- The Header contains a description of the input file. -->
 <Author> Peter Greck </Author>                <!-- Add here the name of the author of this input file. -->
 <Version> 1.0 </Version>                      <!-- Add here the version number      of this input file. -->
 <Content> Terahertz quantum cascade lasers operating up to ~200 K with optimized oscillator strength and improved injection tunneling
           S. Fathololoumi et al.              <!-- Add here some information that describes the content of the input file. -->
           Optics Express 20, 3866 (2012)
           This input file calculates ...
           CPU time: ...

<!-- Here comes the main part, see below... -->



 <Directory                    Comment="Name of output folder for calculated results."         > Output/my_input_file </Directory>

 <WriteOutputEveryNthIteration Comment="Determines how frequent output is written."            > 2                    </WriteOutputEveryNthIteration>

 <MaxNumberOfEigenstates       Comment="Determines number of eigenstates written in output."   > 15                   </MaxNumberOfEigenstates>

If <MaxNumberOfEigenstates> ist not present, 25% of all possible eigenstates are written out. The total number of possible eigenstates corresponds to the total number of spatial grip points. More information on output of eigenstates.

<FormatAsciiEnabled           Comment="Enable output for ASCII files containing data in columns"              > true </FormatAsciiEnabled>
<FormatAsciiExt               Comment="Specify the file extension for ASCII files containing data in columns."> .dat </FormatAsciiExt>

<FormatGnuPlotEnabled         Comment="Enable GnuPlot output based on text output."           > true                 </FormatGnuPlotEnabled>
<FormatGnuPlotExt             Comment="Specify the file extension for GnuPlot output."        > .gnu.plt             </FormatGnuPlotExt>

<FormatVTKPlotEnabled         Comment="Enable VTK output."                                    > true                 </FormatVTKPlotEnabled>

<FormatAvsEnabled             Comment="Enable AVS output."                                    > true                 </FormatAvsEnabled>
<FormatAvsExt                 Comment="Specify the file extension for AVS files."             > .avs                 </FormatAvsExt>
<FormatAvsBinary              Comment="Specify AVS output mode [binary|ascii]"                > binary               </FormatAvsBinary>

Currently, there is no text output included in the code.

<FormatTextEnabled            Comment="Enable text output."                                   > true                 </FormatTextEnabled>  (Currently, there is no text output included in the code.)
<FormatTextExt                Comment="Specify the file extension for text files."            > .dat                 </FormatTextExt>      (Currently, there is no text output included in the code.)


Here, one can define variables (Constant) that are used further below in the input file, e.g. Doping = 1e16. It is useful to add a meaningful comment as well as units for the variables. (Check: Can we also use strings as variables?)

 <Name  Comment="Doping concentration"> Doping </Name>
 <Value Unit="1/cm^3"> 1e16 </Value>

 <Name  Comment="Al concentration of barriers."> Barrier </Name>
 <Value Unit="[0..1]"> 0.15 </Value>
Iterator Example 1

An Iterator can be used to sweep over variables, e.g to do calculations for several temperatures. Here, the iterator accepts a certain list of values.

In this example, a calculation for the temperatures Temperature = 100 [K], 150 [K], ..., 200 [K] is performed.

 <Name    Comment="Temperatures to sweep.">Temperature</Name>
 <Values  Unit="K"> 100, 150, 175, 190, 200 </Values>
Iterator Example 2

The Iterator also accepts a range of values from an Initial value to a Final value with a certain number of Steps. Here, the Final value and the number of Steps is given, and the step width, i.e. the Delta, is derived.

In this example, a calculation from Temperature = 100 [K] to Temperature = 200 [K] is performed for 11 different temperatures.

 <Name    Comment="Temperatures to sweep.">Temperature</Name>
 <Initial Unit="K"> 100 </Initial>
 <Final   Unit="K"> 200 </Final>
 <Steps   Comment="Note that the value is calculated via Initial+[0..Steps-1]*(Final-Initial)/(Steps-1)"> 11 </Steps>
Iterator Example 3

The iterator also accepts a range of values starting from an Initial value and performing a certain number of Steps with a step width of Delta. Here, the step width, i.e. the Delta, and the number of sweep Steps is given, and from these, the final value is derived.

In this example, a calculation for 11 different temperatures with a step width of 10 [K], starting from Temperature = 100 [K] is performed.

 <Name    Comment="Temperatures to sweep.">Temperature</Name>
 <Initial Unit="K"> 100 </Initial>
 <Delta   Unit="K"> 10 </Delta>
 <Steps   Comment="Note that the value is calculated via Initial+[0..Steps-1]*Delta">11</Steps>

Device definition

Temperature as value


 <Temperature Unit="K" > 200         </Temperature>  <!-- Here, a temperature of 200 K is specified. This is the lattice temperature, not the temperature of the electrons. -->

or temperature as a variable.

<Temperature Unit="K" > Temperature </Temperature>  <!-- Here, instead of a number, the variable Temperature is used which has to be specified in the <Variables> section. -->
Specify grid

This is the grid spacing. It determines the number of grid points in the device. The more grid points, the longer the CPU time. If the grid spacing is modified, then always a slightly different structure is calculated as the barrier height and widths are adjusted to match the grid spacing, see <AdjustBandedge>. Even if the grid spacing exactly matches the layer widths for all variations of grid spacings, in all cases slightly different structures are calculated as the dispersion relations change with grid spacing. If the grid is too fine, then one cannot calculate any more sufficiently large energies. In this case, the results would not be reliable.

 <Spacing  Unit="nm"> 0.25 </Spacing>

Orientation and strain

In order to calculate strain related properties, information on the substrate material is required. Strain can be included (yes) or excluded (no). The crystal structure can be Zincblende or Wurtzite. The lattice constants of the substrate are used to calculate strain. The (hkl) Miller indices are needed in order to define the orientation of the substrate on which the heterostructure is grown, i.e. the crystal coordinate system has to be rotated into the simulation coordinate system. The growth direction (simulation axis z) is perpendicular to the substrate (xy plane). The crystal and thus all anisotropic material properties are rotated accordingly.

<Crystal>Zincblende</Crystal>  <!-- Zincblende or Wurtzite -->

 <Material Base="GaAs"/>       <!-- substrate material -->

 <z_axis>    <!-- The heterostructure growth direction is along z (simulation direction). -->
  <h>0</h>   <!-- (hkl) are the Miller indices of the plane perpendicular to the z direction. -->

  <h>0</h>   <!-- (hkl) are the Miller indices of the plane perpendicular to the y direction. -->

<Strain Calculate="no"/>       <!-- include strain (yes/no) -->

Energy grid

Here, the energy grid spacing is defined. The energy grid is homogeneous. The Nodes are the number of energy grid points. The more energy grid points, the longer the CPU time. Example: If one has defined an energy Range of 0.3 [eV], then by choosing 601 Nodes, the resulting energy grid spacing is 0.5 [meV]. If the barrier height in the conduction band is 0.5 [eV], then Range should also be around 0.5. The larger the barrier, the higher the Range.

  <Nodes            Comment="Number of energy grid points."> 601 </Nodes>
  <Range  Unit="eV" Comment="Offset is based on conduction band edge, i.e. Input/BandEdge_conduction_input.dat."> 0.3 </Range>

Definition of layers

<!-- Begin Layers -->
<Width    Unit="nm"> 4 </Width>
<Material Base="GaAs"/>

or alternatively use conduction band edge as defined by <ConductionBandOffset>.

<Material Base="GaAs" CalculateBandedge="no"/>

or (default) calculate conduction band edge by using valence band offset (<ValenceBandOffset>) and temperature dependent band gap, i.e. \(E_\text{c} = E_\text{v} + E_\text{gap}(T)\), taking into account the Varshni parameters.

<Material Base="GaAs" CalculateBandedge="yes"/>
<Probes   Comment="Specifies the relative scattering strength within this layer. Set to zero to remove all probes from this layer.">1.0</Probes>
<Doping   Unit="1/cm^3"> LeadDoping </Doping> <!-- Here, instead of a number, the variable LeadDoping is used which has to be specified in the <Variables> section. -->

It is very useful to define one (or several) periods of a quantum cascade laser with such a marker (BeginCluster). This allows one to define a bias per period very conveniently, i.e. a defined electrostatic potential drop per period. Here, a 4.3 [nm] wide AlGaAs (Al(x)Ga(1-x)As) layer with an alloy concentration of x = 0.15 is defined. The layer is undoped. Here, two labels, namely Center and QCL, are assigned to a cluster which begins here and ends at <Layer EndCluster="Center, QCL">.

<!-- Begin Period -->
<Layer BeginCluster="Center, QCL">
   <Width    Unit="nm"> 4.3 </Width>
<!-- <Material Base="AlGaAs" AlloyX="Barrier" /> <!-- Here, instead of a number, the variable Barrier is used which has to be specified in the <Variables> section. --> -->
     <Material Base="AlGaAs" AlloyX="0.15" />
   <Probes> 1.0 </Probes>
   <Doping   Unit="1/cm^3"> 0 </Doping>

<Layer>                            <!-- Here, a 7.9 nm wide GaAs layer is defined. The layer has a doping concentration of 6e16 1/cm^3. -->
   <Width    Unit="nm"> 7.9 </Width>
   <Material Base="GaAs"/>
   <Probes> 1.0 </Probes>
   <Doping   Unit="1/cm^3"> 6e16 </Doping>

<Layer EndCluster="Center, QCL">   <!-- Here, a 5.5 nm wide GaAs layer is defined. The layer is undoped. -->
   <Width    Unit="nm"> 5.5 </Width>
   <Material Base="GaAs"/>
   <Probes> 1.0 </Probes>
   <Doping   Unit="1/cm^3"> 0 </Doping>
<!-- End Period -->

   <Width    Unit="nm"> 4.0 </Width>   <!-- Here, a 4.0 nm wide GaAs layer is defined. The layer has a doping concentration determined by the variable LeadDoping which has to be specified in the <Variables> section. -->
   <Material Base="GaAs"/>
   <Probes> 1.0 </Probes>
   <Doping   Unit="1/cm^3"> LeadDoping </Doping>
<!-- End Layers -->

Each <Layer> can have an additional flag:

<AdjustBandedge Comment="yes or no">yes</AdjustBandedge>  <!-- default is yes -->

The ratio of the desired barrier width (e.g. 1.1 [nm]) to the width used in the simulation (e.g. grid spacing 0.5 [nm] ==> 2*0.5 [nm] = 1.0 [nm]) is added to the barrier height to keep the area of the barrier the same. This approach compensates the loss in accuracy when using a larger grid spacing very well. This applies analogously to the well width. The effect of this flag can be seen in these files in case they are plotted on top of each other:

  • Input\BandEdge_conduction_input.dat

  • Input\BandEdge_conduction_adjusted.dat


<!-- Begin Contacts -->
  <Lead Name="Source">        <!-- Define a source contact with voltage V = 0.000 V. -->
         <Initial Unit="V"> 0.000 </Initial>
   <Temperature Unit="K"> 300 </Temperature> <!-- Define the temperature of the contact. If not defined, then the device temperature is used -->

  <Lead Name="Drain">        <!-- Define a drain contact with voltage V that is varied during the calculation starting from V = 0e-3 V. -->

(Check: Probably the terms ``Source`` and ``Drain`` are required. This should be checked.)

In the defined structure, there is a region (here the cluster is called Center) where the given voltage drops, e.g. <Cluster>Center</Cluster>. One can include the leads into this region, then the whole voltage drops over the whole structure. Typically, for QCL simulations, the given bias voltage refers to one period only, i.e. the bias voltage per period. In this case, the lead regions are not included. Therefore, the actual voltage drop between the left and the right boundary grid points of the whole structure is larger than the specified value of the voltage drop because the structure is larger than the QCL period as the leads have to be taken into account when calculating the length of the structure.

Example: If one applies 50 mV, this difference corresponds to an energy difference of 50 meV of the conduction band edge of the leftmost and rightmost grid point of the structure, or of the leftmost and rightmost grid point of the region where the bias drops. The effective electric field can be calculated as follows:

field = applied bias / relevant region

Depending on the definition in the input file, the relevant region can include the lead regions or not.

Option 1

A voltage is specified (which typically corresponds to a bias/period).

<Initial Unit="V"> 0e-3  </Initial>
<!-- <Delta   Unit="V"> 3e-3  </Delta> -->  <!-- For the meaning of Initial, Delta, Final and Steps see the documentation of <Iterator>. -->
<Final   Unit="V"> 60e-3 </Final>
<Steps   Unit="" > 15    </Steps>      <!-- If Steps = 1, only one voltage is calculated. If Steps = 0, then automatically, Steps = 1 is assumed. -->
Option 2

As an alternative, instead of specifying a voltage, one can also specify an electric field. Internally, the relevant bias is calculated from the field, and then the calculation starts.

<SpecifiedByElectricField Comment="Specifies the voltage as electric field across the device.
                                   The units of Initial, Final, and Delta are [kV/cm].
                                   Valid input is [yes|no]. Default is [no]."> no </SpecifiedByElectricField>
<Initial Unit="kV/cm"> 0  </Initial>
<!-- <Delta   Unit="kV/cm"> 5  </Delta> -->  <!-- For the meaning of Initial, Delta, Final and Steps see the documentation of <Iterator>. -->
<Final   Unit="kV/cm"> 30 </Final>
<Steps   Unit=""     > 5  </Steps>


<!-- End Contacts -->

MSB model definition - Multi-scattering Büttiker probe model

MSB single-band (1Band = single-band)


Piezoelectric and pyroelectric charge densities can be included (yes) or excluded (no) in the Poisson equation.

<Poisson Piezo="no" Pyro="no"/>

Define numerical parameters

<MaxGreenIts                             Comment="Max. iterations within the G^R cycle."> 25 </MaxGreenIts>
<MaxGreenDelta            Unit="1/nm/eV" Comment="Max. change between two cycles to abort iteration."> 3e-7 </MaxGreenDelta>

<MaxProbeIts                             Comment="Max. iterations for current conservation calculation. Only used if Core.ProbeMode=iterative"> 15 </MaxProbeIts>
<MaxProbeNorm             Unit="A/cm^2"  Comment="Max. absolute leakage current to abort iteration."> 1e-9 </MaxProbeNorm>

<MaxPoissonOuterIts                      Comment="Max. outer Poisson iterations where G^R is recalculated."> 25 </MaxPoissonOuterIts>

<MaxPoissonOuterIts> is an important number. Make sure that this number is not exceeded. If it is exceeded, the calculation did not converge! In the .log file, the following information is written to inform about the number of Poisson iterations: Poisson iteration 30

<MaxPoissonInnerIts                      Comment="Max. inner Poisson iterations where the density predictor is used."> 15 </MaxPoissonInnerIts>

<MaxPoissonDensityNorm    Unit="1/cm^3"  Comment="Max. absolute density   deviation to abort Poisson iterations."> 1e9 </MaxPoissonDensityNorm>
<MaxPoissonCorrectorNorm  Unit="V"       Comment="Max. absolute potential deviation to abort Poisson iterations."> 1e-4 </MaxPoissonCorrectorNorm>

Now we define the core properties of the MSB method.

 <BallisticCalculation Comment="yes or no">no</BallisticCalculation>  <!-- default is no -->

In a ballistic calculation, where scattering is not taken into account, the following parameters are set to zero:

  • Probes

  • ScatteringStrengthConst

  • ScatteringStrengthBP

  • ScatteringStrengthMSB

  • bulk (default): The initial guess of the electrostatic potential assumes a bulk semiconductor where the resulting density is neutral at the grid point 0 or n, respectively.

  • zero: The initial guess of the electrostatic potential is zero.

    <ProbeMode   Comment="Specify method to calculate current conservation. Novel (more stable) method is 'direct'. [direct|iterative]"> direct </ProbeMode>

Note that direct and iterative can lead to different results. They are not equivalent.

  • The iterative model is explained in Section 5.2 in the PhD thesis of Peter Greck. Here, for each probe, virtual chemical potentials \(\mu_\text{p}(z)\) are calculated. The units are [eV]. Therefore, for each probe a virtual chemical potential exists which leads to a Fermi distribution of each probe \(f(E,\mu_\text{p}(z))\). This virtual chemical potential \(\mu\) refers to the total probe, i.e. AC + LO.

  • The direct model is explained in Section 7.2 in the PhD thesis of Peter Greck. Here one assumes that the distribution function \(f_\text{p}(z)\) of each probe is a linear combination of the source and drain distributions, \(f_\text{S}\) and \(f_\text{D}\), respectively, where S means Source, and D means Drain.

    \(f_\text{p}(E) = c_\text{p} f_\text{S}(E) + (1 - c_\text{p}) f_\text{D}(E)\)

    The coefficients \(c_\text{p}\) for each probe are in the interval [0,1] and are dimensionless. Here, a linear system of equations has to be solved to obtain the coefficients \(c_\text{p}\). In contrast to the iterative model, the virtual chemical potentials \(\mu_\text{p}\) for the probes cannot be extracted in this case. For the results presented in the PhD thesis of P. Greck, always the direct model was used.

    <VoltageMode Comment="Specify handling of applied bias voltage. Drop-mode ensures that the specified voltage drops along the device. Flat-mode uses Neumann boundary conditions that do not guarantee a complete voltage drop. [drop|flat]"> drop </VoltageMode>

The total scattering strength that is used for the calculation is the sum of the following three products, i.e. total = Model #1 + Model #2 + Model #3:

<Probes> * <ScatteringStrengthMSB> +

<Probes> * <ScatteringStrengthBP> + (This term is usually zero.) ==> These are the “normal” Büttiker probes (see [DattaETMS1995]).

<Probes> * <ScatteringStrengthConst> +  (This term is usually zero.) ==> These are the “oldest” type of Büttiker probes (see [DattaETMS1995]).

We recommed to only define <ScatteringStrengthMSB> and set <ScatteringStrengthBP> and <ScatteringStrengthConst> to zero.

Model #1 for scattering

Here, we can define a scattering strength for both LO, and AC phonon scattering.

<ScatteringStrengthMSB      Unit="[0..1]" Comment="Novel scattering via MSB. Scattering strength is calculated from material parameters. This parameter tunes the scattering rates from 0.0=disabled over 1.0=normal to X=amplified."> 1.0 </ScatteringStrengthMSB>
<ScatteringStrengthMSB_AC   Unit="[0..1]" > 1.0 </ScatteringStrengthMSB_AC>
<ScatteringStrengthMSB_LO   Unit="[0..1]" > 1.0 </ScatteringStrengthMSB_LO>
Model #2 for scattering

These are the normal Büttiker probes (see [DattaETMS1995]).

<ScatteringStrengthConst Unit="eV"     Comment="Additional constant scattering strength for every probe. Typical values are 0.001-0.01 and zero disables this scattering mechanism."> 0.0 </ScatteringStrengthConst>
Model #3 for scattering

These are the oldest type of Büttiker probes (see [DattaETMS1995]).

<ScatteringStrengthBP    Unit="eV"     Comment="Additional scattering via Büttiker-Probes. Typical values are 0.001-0.01."> 0.0 </ScatteringStrengthBP>



  <Cluster                       Comment="Specfiy the cluster where optical gain should be calculated."> Center </Cluster>
  <PhotonEnergyInitial Unit="eV" Comment="Min. photon energy for gain calculation."                    >  3e-3  </PhotonEnergyInitial>
  <PhotonEnergyFinal   Unit="eV" Comment="Max. photon energy for gain claculation."                    > 30e-3  </PhotonEnergyFinal>
  <PhotonEnergySteps   Unit=""   Comment="Number of photon energies for gain calculation."             > 100    </PhotonEnergySteps>


The more PhotonEnergySteps are used, the more ragged the gain curve gets. Decreasing this number, smoother gain curves can be obtained. This is normal. If the photon energy grid spacing is too small, one tries to resolve energies that are not resolved within the calculated Green’s functions that are used in the gain calcuations, i.e. one tries to calculate several photon energies within one discrete energy step of the Green’s function. An obvious solution would be to increase the energy grid resolution of the Green’s functions. However, it is recommended to rather prefer a course photon energy grid as a high photon energy grid resolution does not lead to significantly more insight.


  • One can disable the scattering within the barriers. This can reduce computational time with very minor influence on overall results.

     <Probes> 0.0 </Probes>  <!-- disabled -->
     <Probes> 1.0 </Probes>  <!-- enabled  -->

    Important: If there are no probes inside the barrier, the iterative calculation of the probes is much more stable while hardly changing the physics because the occupation in the barriers is basically zero. Moreover, this avoids the numerical explosion of the linear system of equations through a division by zero which leads to NaN (not a number) in the .log file.

  • How can I do a ballistic calculation? There are three ways to achieve this.

    1. Use

      <BallisticCalculation Comment="yes or no"> yes </BallisticCalculation>
    2. Set all <Probes> to zero, i.e. 0.0.

    3. Set all of the following parameters to zero: <ScatteringStrengthMSB>, <ScatteringStrengthBP> and <ScatteringStrengthConst>, i.e. 0.0.

  • For a quantum cascade laser simulation it can make sense to simulate 5 or 7 periods and take the middle period as the one where one extracts the results like local density of states, electron density, current density, gain, …

  • If one simulates only one period of a QCL, then one should add a doping layer at the left and right boundary (e.g. of width 4.5 nm). This should not be necessary if 5 periods or more are simulated. In this case, doping could be omitted.