In this example we will see how to perform parallel calculations to study electronic and optical properties of a InGaN/GaN Quantum Well (QW) structure with random alloy InGaN active region, using Empirical Tight Binding (ETB).

With tibercad (linux version), it is now possible to run a simulation with a MPI parallel execution on a multicore machine with the command:

tibercad -n number_mpi_proc input_file.tib

where number_mpi_proc is the number of processes to run in parallel.
For example, on a 4-core machine, the command

tibercad -n 4 input_file.tib

will parallelize the simulation on 4 cores of the local processor.

When running in parallel with MPI, the default is to consider one single device and run solvers in parallel for this device.
Here instead we show, as an alternative, how to run in parallel a simulation on a given number of different instances of the device, with or without parallel solvers.
This can be useful when considering statistical sets of random alloy configurations.

In Device section, we can set some options controlling parallelization inside the Parallel block:

Parallel
 {

 # one device (random alloy sample) for each process
 mpi_processes_per_device = 1 # 2

 # let's do FEM calculations in serial
 mpi_processes_per_mesh = 1

 }

Here, mpi_processes_per_device sets the number of processes assigned to a given device instance. By default mpi_processes_per_device  is set to the value of  number_mpi_proc, that is a single device is parallelized on the chosen number of processes.

If mpi_processes_per_device < number_mpi_proc, then more than one random alloy device is calculated (in serial or parallel).

Thus, for example, on a 4-core machine, with:

tibercad -n 4 file.tib 
 mpi_processes_per_device = 2

2 devices are calculated, each with parallelization on 2 cores.

With

tibercad -n 4 file.tib 
 mpi_processes_per_device = 1

then 4 devices are calculated in parallel, one on each of 4 cores.

Finally, by setting:

mpi_processes_per_mesh = 1

we impose to execute FEM calculations in serial in any case.
Atomistic calculations, which are much more time demanding, can be instead parallelized, depending on the case.

NB: In the current version it is not possible to perform a parallel simulation of a device where different portions of the mesh are assigned to different FEM Modules.
Therefore, a simulation can be parallelized only if ALL the FEM Modules are applied to the
same set of Physical Regions

In order to execute correctly this example you should have the following files in the same working directory:
led_TB_random_alloy.tib: input file for TiberCAD
led_TB_random_alloy.msh: mesh file produced by the GMSH script led_TB_random_alloy.geo
dd_bias_2.6.tsv: saved solution for V=2.6V

DEVICE STRUCTURE

The device structure is defined in the geometry .geo file and is the following:

• 150 nm n-doped GaN buffer
• 3 nm InGaN QW with 25% In content
• 4 nm undoped GaN barriers
• 30 nm AlGan ebl
• 50 nm p-doped GaN cup and 50 nm p+ GaN buffer

In the Device section, the QW heterostructure is described: the InGaN quantum well region, the GaN barriers and buffer regions.

The crystal directions are defined as following for a wurtzite crystal structure:

y-growth-direction = (1,0,-1,0)

z-growth-direction = (-1,2,-1,0)

x-growth-direction = (0,0,0,1)

The regions are defined in the usual way, then a cluster named quantum for k.p calculations and a cluster atomistic for ETB calculations are defined, in this way:

# for k-dot-p
 Cluster quantum
 {
 regions = (well1, qbarrierl, qbarriert, buffer_q)
 }

 # for tight binding
 Cluster atomistic
 {
 regions = (well1, qbarrierl, qbarriert)
 }

ATOMISTIC STRUCTURE

An atomistic representation of the above defined atomistic cluster is generated by means of the Atomistic block

Atomistic tb
 {

 reference_region = qbarrierl
 regions = atomistic
 passivation = yes
 print = (xyz)

 random_alloy = true

 random_generator_seed = MPI_DEV_KEY # 1

 supercell_size_z = 20 # size in Angstrom
 supercell_size_y = 20

}

With

random_alloy = true

we set a random configuration of the InGaN alloy, instead of the default VCA approach.
In this way, the atoms of the components species are randomly distributed in the alloy region and their parameters are not averaged.

reference_region = qbarrierl

The reference region is chosen to provide the lattice parameters with which the crystalline structure is built. In this case the lattice is that of GaN, the material composing the qbarrier.
In the QW region, the In atoms are then substituted in the lattice basis.

By default a 2D periodicity is applied in yz-plane orthogonal to the x growth direction of this 1D QW structure.

In this case, however, we define a custom supercell, with a 20x20 Angstrom size.

supercell_size_y = 20

supercell_size_z = 20

Therefore, the appropriate periodicity vectors will be applied to this defined supercell, not to the default minimal cell.

Passivation is finally performed at the ends of the heterostructure

passivation = yes

A print instruction gives in output the atomic structure for a visualization.

Here we add also an optional seed for the random generator, in order to allow for repeatability and reload of results

random_generator_seed = MPI_DEV_KEY  # 1

In this way, if in the future one wants to repeat or extend calculations with the same random alloy configuration, one has just to set random_generator_seed to the value used in previous calculations, e.g.

random_generator_seed = 1

The variable MPI_DEV_KEY  is useful when more than one device is calculated in parallel. In this case, each device (that is its random alloy configuration) can be identified by the value of MPI_DEV_KEY  (see RUN SIMULATION below).

SIMULATION MODULES

We first perform elasticity calculations to obtain the strain tensor in the heterostructure and to apply a strain deformation to the defined atomistic structure.
This result provides a reasonable first guess for the valence force field (VFF) solver, which will yield the final relaxed structure.
We then perform the drift-diffusion model to get the bias point solution for the desired potentials and band profiles.
Here we can also load a pre-calculated solution.
Quantum ETB calculations are finally performed on the VFF-relaxed structure, to get the electron and hole states in the QW and the optical recombination spectrum.

1. Elasticity

Note that it is always highly advisable to solve first elasticity, before a VFF relaxation. This will provide a reasonable first guess for the VFF solver and also will apply the correct deformation to the mesh, assuring that all the atoms are still associated to the correct FEM elements.

We solve elasticity for the lattice-mismatch induced strain

Module elasticity 
{

 name = strain
 regions = all

 plot = (Strain, Stress)

 # deform mesh and atomistic structure
 # this is needed here since the structure is built with GaN
 # as reference crystal
 
 mesh_deformation = true


 Physics 
 {
 body_force lattice_mismatch 
 {
 reference_material = GaN
 y-growth-direction = (1,0,-1,0)
 z-growth-direction = (-1,2,-1,0)
 x-growth-direction = (0,0,0,1)
 }
 
 }

 Contact substrate  
 {type = clamp}

}

Note that, in order to apply a strain deformation to the defined atomistic structure, it is convenient to define a mesh deformation with the default number of shape_iterations.

mesh_deformation = true

In this way, strain will be computed iteratively until the convergence on the structure deformation is reached.

2. VFF

Now, based on the displacements obtained from the elasticity solution and projected onto the atomistic structure, we proceed to its relaxation with a valence force field (VFF) approach. We thus solve VFF on the atomic structure:

Module vff
{

 regions = atomistic 
 atomistic_structure = tb

 # fixes first layer of atoms on x-plane,
 # with a certain tolerance on x coordinate
 boundary_conditions = substrate
 substrate_plane = x
 substrate_tol = 0.519

 plot = StrainCells

 Physics
 {

 keating automatic
 {
 use_four_wz_parameters = true
 }
 }
}

with the b.c.

boundary_conditions = substrate

that is, first layer of atoms on x-plane is fixed.

3. Drift-diffusion

As for drift-diffusion, as usual we define a simulation

name = driftdiffusion

belonging to the model driftdiffusion and associated to the whole device (deafult choice).

In this example, we suppose to have performed previously a sweep to the bias point V=2.6.
So, here we just reload the solution state at the bias point of interest

load_state = "dd_bias_"BIAS".tsv"

4. Empirical Tight-Binding

For the atomistic quantum calculations, we define an empirical_tb simulation, named tb:

Module empirical_tb
{
 regions = atomistic
 name = tb

 atomistic_structure = tb
 Harrison_scaling = true # set to false if no strain is calculated

 potential_simulation = dd
 
 # plot = (tbstates, MeshStates)

 Solver
 {
 
 # use dynamic_search to improve the solver first guess 
 dynamic_search = true 

 # dynamic_offset = 0.1 # default 0.1

 num_conduction_eigenvalues = 8
 num_valence_eigenvalues = 8 

 long_tolerance = 1e-4

 }

}

In the Module empirical_tb we define the associated regions, given by

regions = atomistic

ETB Hamiltonian will be created and solved based on the atomistic structure named tb which has been previously built in the Atomistic section and relaxed through the elasticity model and the valence force field one.

atomistic_structure = tb

Information about potential profile to be applied (due to built-in, polarization fields, etc.) is given by

potential_simulation = dd

In the Solver block, we define the number of eigenstates to be calculated in valence and conduction band:

num_valence_eigenvalues = 8

num_conduction_eigenvalues = 8

Moreover, here we set to true the option to use dynamic_search, to improve the solver first guess:

dynamic_search = true

ETB calculations will be applied to a relaxed atomistic structure, where the atomic positions are displaced from the reference ones. Since we want the TB parameters to be scaled with respect to this displacement, according to Harrison law, we define

Harrison_scaling = true

In this way, the effect of strain relaxation in the heterostructure will be correctly talken into account in ETB calculations.

4. Optical properties with Empirical Tight-Binding

For the calculation of optical properties, we define a opticstb simulation:

Module opticstb
{
 name = opticstb
 regions = atomistic 

 plot = (optical_spectrum, matrix_elements)

 # put it to write output in 1D
 output_format=grace

 initial_state_model = tb
 final_state_model = tb
 skip_solve = true # skip tb model, use tb results from previous step

 Emin = 2
 Emax = 4 
 dE = 0.001
 broadening = 0.01
 line_shape = gaussian

 k_integration
 {
 k_max = 1 
 number_of_elements = (1,1)
 quadrature_type = trapezoidal 

 refine_k_space = false 
 }

}

The optical matrix and the optical spontaneous emission spectrum are calculated based on ETB Hamiltonian.

Simulation is performed with an integration in k-space, choosing the keyword optical_spectrum in the plot command:

plot = (optical_spectrum, matrix_elements)

k-space integration is defined by the block:

k_integration
 {
 k_max = 1 
 number_of_elements = (1,1)
 quadrature_type = trapezoidal 

 refine_k_space = false 
 }

RUN SIMULATIONS

We may now run tiberCAD to calculate and apply strain deformation first with a macroscopic elasticity model (strain) and then with a valence force field one (vff), we solve then driftdiffusion (dd) to get a bias point solution and ETB for calculation of eigenvalues of holes and electrons (tb). Finally, optical spectrum is obtained through opticstb.

solve = (strain, vff, dd, tb, opticstb)

Usually one may choose to perform a set of simulations of this kind, to get a statistical ensemble of random structures.
This can be done e.g. in parallel.

In this example, since we have set, in Parallel block:

 mpi_processes_per_device = 1

by executing tibercad with the command:

tibercad -n 4 file.tib 

4 different devices are simulated in parallel, each one with a different random alloy choice.

In this case the variable MPI_DEV_KEY is used as an identifier for a single device process group:

resultpath = "output_tb_"MPI_DEV_KEY

This means that in general the output will consist of several output directories (in this case 4), one for each "device" calculated.
Each one will show eigenvalues and eigenstates slightly different, due to the different configuration of In atoms in the InGaN random alloy.

OUTPUT

By using jmol, (an open-source Java viewer for chemical structures in 3D, see http://jmol.sourceforge.net) we can load the atomic structure generated by tiberCAD Atomistic Generator, contained in the file tb.xyz (one in each output directory).
It is made by a GaN/InGaN/GaN heterostucture grown along the x-axis direction. Note that a supercell with 20 Angstrom size in yz plane has been created. The supercell structure is then periodical in the yz plane. A different atomistic structure can be seen for each device instance, due to the different random alloy configuration.

As for tight-binding simulation, the file tb.dat contains a table with the calculated eigenvalues for electrons and holes, together with their occupation index. Finally, opticstb.dat contains the calculated optical matrix elements and the file opticstb_spectrum.dat contains the spontaneous power densities Px, Py and Pz.

ATTACHMENTS

• led_TB_random_alloy.geo
  
• led_TB_random_alloy.msh
  
• led_TB_random_alloy.tib
  
• dd_bias_2.6.tsv