5.3. Input file syntax

Input file

The structure to be simulated as well as all simulation parameters have to be specified in a text file by the user. In addition to the XML format, the nextnano++-style format with the extension .negf is supported (since 2022 Apr).

Both formats of the sample input files are provided in the installation folder, e.g. THz_QCL_Fathololoumi_OpticsExpress2012.xml and THz_QCL_Fathololoumi_OpticsExpress2012.negf. This documentation explains the syntax in the .negf format, but the keywords and hierarchical structure are identical between two formats except for <Variables> (difference described below).

You can convert the existing XML input files to the new format using nextnano.MSB executables later than 2022 Apr. If additional command line option -p (or --parse) is set, the execution stops after the input file and database have been read and parsed. In case XML files were read in, tentative conversions into the new format are output.

General syntax

nextnano.MSB{         # Version = 1.1.0   # This number must be consistent to the version of the executable used.

In the XML format,

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

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

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

Comment Section

Header{                                     # The Header contains a description of the input file.
   Author = "Peter Greck"                   # name of the author of this input file
   Version = 1.0                            # version number of this input file
   Content = "See comments in input file."  # information that describes the content of the input file
   #       Terahertz quantum cascade lasers operating up to ~200 K with optimized oscillator strength and improved injection tunneling
   #       S. Fathololoumi et al.
   #       Optics Express 20, 3866 (2012)
   #         ...
   #       CPU time: ...
   }

   # Here comes the main part, see below...
}

Output

Output{
   Directory                    = Output/my_input_file    # Comment = "Name of output folder for calculated results."

   WriteOutputEveryNthIteration = 2                       # Comment = "Determines how frequent output is written."

   MaxNumberOfEigenstates       = 15                      # Comment = "Determines number of eigenstates written in output."

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   = true         # Comment = "Enable output for ASCII files containing data in columns"
FormatAsciiExt       = .dat         # Comment = "Specify the file extension for ASCII files containing data in columns."

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

FormatVTKPlotEnabled = true         # Comment = "Enable VTK output."

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

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

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

Variables

New: In the .negf format, variables must be defined with the dollar sign as in nextnano++ (see Input Syntax and sample input file in the installation folder).

Iterator feature is realized by the nextnanomat feature Template or nextnanopy Sweep. This simply generates multiple input files, which can be run in parallel (see Options: Simulation).

$Temperature = 300   # (DisplayUnit:K) Temperature to sweep. (ListOfValues:100,150,175,190,200) (RangeOfValues:From=100,To=200,Step=11)
$Doping = 1e16       # (DisplayUnit:1/cm^3) Doping concentration
$Barrier = 0.15           # Al concentration of barriers.

Below is the documentation of <Variables> in the old XML format input file.

In the XML format, one can define variables (Constant) that are used further below in the input file, e.g. Doping = 1e16.

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

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

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.

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

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.

<Iterator>
   <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>

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.

   <Iterator>
      <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>
   </Iterator>
</Variables>

Device definition

Temperature as value

Device{

   Temperature = 200            # Unit = K  # This is the lattice temperature, not the temperature of the electrons.

or temperature as a variable.

Temperature = $Temperature   # Unit = K  # Here, instead of a number, the variable $Temperature is used.

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.

Grid{
   Spacing = 0.25     # Unit = nm
}

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   # Zincblende or Wurtzite

Substrate{
   Material{
      Base = GaAs      # substrate material
   }
}

Orientation{
   z_axis{    # The heterostructure growth direction is along z (simulation direction).
      h = 0   # (hkl) are the Miller indices of the plane perpendicular to the z direction.
      k = 0
      l = 1
   }

   y_axis{
      h = 0   # (hkl) are the Miller indices of the plane perpendicular to the y direction.
      k = 0
      l = 1
   }
}

Strain = 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.

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

Definition of layers

# Begin Layers
Layer{
   Thickness = 4   # Unit = nm
   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 = 1.0          # Comment = "Specifies the relative scattering strength within this layer. Set to zero to remove all probes from this layer."
   Doping = $LeadDoping  # Unit = 1/cm^3  # Here, instead of a number, the variable $LeadDoping is used.
}

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"
   Thickness = 4.3           # Unit = nm

   # Material{
   #    Base = AlGaAs
   #    AlloyX = $Barrier    # Here, instead of a number, the variable $Barrier is used.
   # }

   Material{
      Base = AlGaAs
      AlloyX = 0.15
   }
   Probes = 1.0
   Doping = 0        # Unit = 1/cm^3
}

Layer{               # Here, a 7.9 nm wide GaAs layer is defined. The layer has a doping concentration of 6e16 1/cm^3.
   Thickness = 7.9   # Unit = nm
   Material{
      Base = GaAs
   }
   Probes = 1.0
   Doping = 6e16     # Unit = 1/cm^3
}

Layer{               # Here, a 5.5 nm wide GaAs layer is defined. The layer is undoped.
   EndCluster = "Center, QCL"
   Thickness = 5.5   # Unit = nm
   Material{
      Base = GaAs
   }
   Probes = 1.0
   Doping =   0      # Unit = 1/cm^3
}
# End Period

Layer{
   Thickness = 4.0   # Unit = nm   # Here, a 4.0 nm wide GaAs layer is defined. The layer has a doping concentration determined by the variable $LeadDoping.
   Material{
      Base = GaAs
   }
   Probes = 1.0
   Doping = $LeadDoping    # Unit = 1/cm^3
}
# End Layers

Each Layer{ can have an additional flag:

Layer{
   AdjustBandedge = yes   # Comment = "yes or no"  # 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

Contacts (leads)

# Begin Lead
Lead{
   Name = Source       # Define a source contact with voltage V = 0.000 V.
   Voltage{
      Initial = 0.000  # Unit = V
   }
   Temperature = 300   # Unit = K # 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.
   Voltage{
      ...
   }

(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).

Voltage{
   Cluster = Center
   Initial = 0e-3      # Unit = V
   # Delta = 3e-3      # Unit = V   # For the meaning of Initial, Delta, Final and Steps see the documentation of <Iterator>.
   Final   = 60e-3     # Unit = V
   Steps   = 15        # Unit = ""  # 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.

      Voltage{
         SpecifiedByElectricField = yes   # 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 = 0   # Unit = kV/cm
         # Delta = 5   # Unit = kV/cm     # For the meaning of Initial, Delta, Final and Steps see the documentation of <Iterator>.
         Final   = 30  # Unit = kV/cm
         Steps   = 5   #  Unit = ""
      }
   } # End Lead
} # End Device

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

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

MSB_1Band{

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

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

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

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      = 15                   # Comment = "Max. inner Poisson iterations where the density predictor is used."

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

Now we define the core properties of the MSB method.

Core{
   BallisticCalculation = no   # Comment = "yes or no"  # 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

  PotentialInit = bulk
  

• 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 = direct  # Comment = "Specify method to calculate current conservation. Novel (more stable) method is 'direct'. [direct|iterative]"
  

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 = drop   # 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]"
  

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    = 1.0  # 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."
ScatteringStrengthMSB_AC = 1.0
ScatteringStrengthMSB_LO = 1.0

Model #2 for scattering

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

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

Model #3 for scattering

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

   ScatteringStrengthBP    = 0.0   # Unit = eV    # Comment = "Additional scattering via Büttiker-Probes. Typical values are 0.001-0.01."
} # End Core

Gain

   Gain{
      Cluster             = Center               # Comment = "Specfiy the cluster where optical gain should be calculated."
      PhotonEnergyInitial = 3e-3    # Unit = eV  # Comment = "Min. photon energy for gain calculation."
      PhotonEnergyFinal   = 30e-3   # Unit = eV  # Comment = "Max. photon energy for gain claculation."
      PhotonEnergySteps   = 100     # Unit = ""  # Comment = "Number of photon energies for gain calculation."
   }
} # End MSB_1Band

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.

Hints

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

  Layer{
     Probes = 0.0   # disabled
  }
  
  Layer{
     Probes = 1.0   # 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

  Core{
     BallisticCalculation = yes   # Comment = "yes or no"
  }
  

  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.