Silicon p-n junction

In this tutorial you will:

  • Learn how to build a p-n junction device.
  • Learn how to study the current-voltage characteristic of such a device.
  • Compare two different methods: ATK-SE and ATK-DFT.
  • Analyze the electronic structure of the p-n junction by studying and plotting the device density of states at zero bias and at reverse and forward bias.

Silicon bulk: Slater-Koster vs DFT-MGGA

  1. Open VNL and create a new project by clicking on Create New. Give the project a Title (here: “Silicon”), select a folder for the project, click OK and then Open to start the project.

    ../../_images/110.png ../../_images/24.png
  2. In the VNL main window, click on the builder_icon icon to open the Builder.

  3. In the Builder, go to the Stash, click Add ‣ From Database and from the menu navigate to Databases ‣ Crystals and search for “silicon” in the search window. Then, select “Silicon (alpha)” and click the plus_icon icon, or double-click its line in the list, to send it to the Stash.



Slater-Koster calculation

You will now set up an ATK semiempirical (ATK-SE) calculation using the Slater-Koster model Hamiltonian.

  1. Send the bulk configuration from the builder_icon Builder to the script_generator_icon Script Generator using the sendto_icon button.

  2. Add the following blocks to the Script by double-clicking the correspoding icons in the Blocks panel:

    • calculator_icon New Calculator.
    • analysis_icon Analysis ‣ Bandstructure.

  3. Open the calculator_icon New Calculator by double-clicking it.

    • Select the ATK-SE Slater-Koster calculator.
    • Use a 25x25x25 k-point sampling.
    • Uncheck the “No SCF iteration” box.
    • Select the “Bassani.Si” basis set.

  4. In analysis_icon Bandstructure, select the “L, G, X, U, K, G” Brillouin zone route and 401 points per segment.

  5. Back in the script_generator_icon Script Generator, set the output filename to, and send the script to the job_manager_icon Job Manager using the sendto_icon button.

  6. Save the script in the window that shows up and use the jm_play_enabled_icon botton to run the job.

  7. Once the calculation is done (it will only take a few seconds) you can find in the LabFloor your Bandstructure object, which is saved in


    Highlight the object and use the Bandstructure Analyzer plugin to plot the bandstructure.


  8. You can zoom into a region close to the band edges and precisely measure the calculated indirect band gap, 1.186 eV.



DFT-MGGA calculation

The TB09 MGGA approximation [aTB09] available in ATK-DFT is fully explained in the Meta-GGA and 2D confined InAs tutorial. In order to compare the Slater-Koster results with results from the MGGA method, you will fit a parameter (c) such the MGGA band gap of silicon matches that calculated with Slater-Koster. Throughout this tutorial you will use this fitted c parameter for the MGGA calculations.

  1. From the builder_icon Builder, send the Silicon bulk configuration to the script_generator_icon Script Generator and add a calculator_icon New Calculator to the script. Then edit the calculator parameters:

    • Select ATK-DFT and the MGGA exchange-correlation functional.
    • Use a 25x25x5 k-point sampling.
    • Select the SingleZetaPolarized basis set in order to speed up calculations.

  2. Add also a analysis_icon Bandstructure analysis object and set it up as you did for the Slater-Koster calculation.

  3. Back in the script_generator_icon Script Generator, set the output filename to

  4. In order to fit the c parameter to the Slater-Koster band gap you will need to run a few calculations for different c parameters and make a linear interpolation. Here, you will run four MGGA calculations for c = 0.9, 1.0, 1.1 and 1.2. A little manual editing of the script is needed:

    • Send the script to the editor_icon Editor using the sendto_icon button, and modify the exchange-correlation section as shown in the image below. Also, add a loop over the c values to run the four calculations with a single Python script (you can download a copy of the script here: script):
  5. Send the script to the job_manager_icon Job Manager and run the job.

  6. On the LabFloor you will now have four different Bandstructure objects included in the “”.

  7. You will use the following script to fit the c-parameter: Save it in the project directory.

  8. Execute the fitting script using the job_manager_icon Job Manager. The plot shown below should pop up.


It is clear that the TB09 c-parameter has a significant influence on the calculated band gap, and that the dependence is quite linear. The script output shows that the fitted value of c is 1.048.


If you take a closer look at the script “”, you will see that band gaps are calculated from a Bandstructure object like this:

bandstructure = nlread("", Bandstructure)[0]
gap_direct   = bandstructure._directBandGap().inUnitsOf(eV)
gap_indirect = bandstructure._indirectBandGap().inUnitsOf(eV)
print "direct gap:   %.2f eV" % gap_direct
print "indirect gap: %.2f eV" % gap_indirect

Silicon device

You will now set up the silicon device.

  1. Open the builder_icon Builder, highlight “Silicon (alpha)” in the Stash, and start creating a Si(100) surface:

    • Navigate to Builders ‣ Surface (Cleave).
    • In the Surface (Cleave) window, choose the [1,0,0] Miller indices and click “Next”.

  2. Keep the 1\(\times\)1 surface unit cell and click “Next”.

  3. Select a thickness of 52 layers and click “Finish”.



    The 52 layers will constitute the central region of your device. You need such a relatively long device (~14 nm) because of the long screening length in the silicon semiconductor, as explained in the Doping in ATK: the Ni silicide - Si interfaces tutorial. Still, as you will see later, this device is not long enough to have perfectly converged results. However, you will use the 52 layers device throughout this tutorial to be able to perform the MGGA calculations faster.

  4. Your Si(100) configuration is now constructed, and you will use it to construct a device with the transport direction along (100). Navigate to Device Tools ‣ Device From Bulk, keep the default parameters, and simply click “OK”.


  5. It will prove convenient to assign tags to atoms in the regions that should be doped p-type and n-type, respectively.

    • Highlight the device configuration in the Stash and go to Selection Tools ‣ By Expression.
    • Type in c < 0.5 and press Enter to select all atoms in the left hand side of the device.
    • Use Selection Tools ‣ Tags to assign the tag “p_doping” to the selected atoms.
    • Follow the same procedure to select atoms with c > 0.5 and assign the tag “n_doping” to those atoms.


Slater-Koster IV curve

You will use the Slater-Koster model to calculate the I-V characteristics for the p-n doped silicon junction.

  1. Send the device configuration from the builder_icon Builder to the script_generator_icon Script Generator.

  2. Add a calculator_icon New Calculator to the script and set the following parameters:

    • Select the ATK-SE Slater-Koster calculator.
    • Use a 7x7x100 k-point sampling.
    • Select the Bassani.Si basis set and uncheck the “No SCF iteration” box.
  3. You can now add all the analysis objects needed. Add the following ones:

    • analysis_icon ElectronDensity
    • analysis_icon DeviceDensityofStates
    • analysis_icon ElectrostaticDifferencePotentials
    • analysis_icon IVCurve
  4. Edit the analysis_icon DeviceDensityofStates parameters:

    • Increase considerably the k-point sampling. Use a 21x21 grid.
    • Energy: from -2 to 2 eV with 401 points.
  5. Edit the analysis_icon IVCurve parameters:

    • Voltage Bias: sample from -1 V to 1 V with 21 points.
    • Energy: from -1.1 eV to 1.1 eV with 401 points.
    • k-point sampling: 21x21.


    A full explanation of all the parameters is given in the Doping in ATK: the Ni silicide - Si interfaces and in the ATK Reference Manual.

  6. Change the output filename to something like

Doping the silicon junction

You need to set up the proper doping to finally have a p-n junction. To this end, send the script to the editor_icon Editor and locate the part of the script that defines the device configuration:

device_configuration = DeviceConfiguration(
[left_electrode, right_electrode]

Just before those lines, add the following lines to define a p-n junction with a doping concentration of 2\(\cdot\)1019 cm-3:

# -------------------------------------------------------------
# Add Doping
# -------------------------------------------------------------
doping_density = 2e+19
# Calculate the volume  and convert it to cm^-3
# Note: right and left electrodes have the same volume,
volume = right_electrode_lattice.unitCellVolume()
volume = float(volume/(0.01*Meter)**3)
# Calculate charge per atom
doping = doping_density * volume / len(right_electrode_elements)
# Add p- and n-type doping to the central region
external_potential = AtomicCompensationCharge([
     ('p_doping', -1*doping), ('n_doping',doping)
# Add doping to left and right electrodes
   AtomicCompensationCharge([(Silicon, -1*doping)]))
   AtomicCompensationCharge([(Silicon, doping)]))

Running the calculation

Send the script to the job_manager_icon Job Manager and run the calculation.

Depending on the analysis objects included in the script, and on the number og parallel processes executed by the Job Manager, the calculation may take a few minutes for the zero bias part and up to few hours to calculate the whole IV curve. The performance of ATK 2014 for different methods (ATK-SE and ATK-DFT) and for a transmission spectrum calculation are reported in the figure below.



The calculations are run in parallel on different Intel Xeon e5472 3.0 GHz machines. The single points at 12 MPI are run with ATK 13.8. Overall, ATK 2014 is 25-30% faster compared to ATK 13.8.

Calculation details: zero-bias, 7x7x100 k-point sampling for SCF part and 21x21 for the transmission spectra. All other parameters are set as in this tutorial.

Note that the Slater-Koster SCF calculation is much faster than the post-SCF calculation of the transmission spectrum.


While the Slater-Koster job is running you can set up the DFT-MGGA calculation. Simply go back to the script_generator_icon Script Generator and change the calculator_icon Calculator block:

  • Choose ATK-DFT.
  • Exchange correlation: MGGA.
  • SingleZetaPolarized basis set.

Remove all analysis objects except the analysis_icon IVCurve. Then set up the device doping and run the calculation.

Analyzing the results

IV Curve

Once the two calculations are done and you have the output nc files in the project directory you will see the two IVCurve objects loaded in the LabFloor:


  1. Select one of the IVCurve objects and use the IV-Plot plugin to plot the IV curve. If you check the Additional plots option, you will also see plots of dI/dV, transmission spectra, and the spectral current.

  2. In order to compare the Slater-Koster and MGGA IV curves, download the script to your project folder and send it to the job_manager_icon Job Manager to create the figure below, where you can compare directly the two IV curves. You can also open the script with the editor_icon Editor and modify it to tune the plot to your best convenience.



From the plot above you can see that the results from the semiempirical Slater-Koster model and the DFT-MGGA method are virtually the same.

Device Density of States

A classical picture of the electronic structure of a p-n junction is the one reported below, where the potential energy is plotted as a function of position along the transport direction of the device.


Fig. 38 Energy band diagram for a p-n junction. Image taken from

Using VNL and ATK it is possible to perform this analysis through the DeviceDensityOfStates or the LocalDensityOfStates analysis objects.


Here, you will analyze the DeviceDensityOfStates and plot it along the device transport direction using a script that you can download here:

In order to run the script and make the plots shown below you need an .nc file containing the analysis objects DeviceConfiguration, DeviceDensityOfStates and ElectrostaticDifferencePotential. You can open the script with the editor_icon Editor and customize the plot according to your wish.

In particular, the script will read:

  • DeviceConfiguration: to determine the spatial resolution along the device by defining a projection list used to plot the DDOS. Also, the electrode voltages are extracted from this object.
  • DeviceDensityOfStates: besides reading and plotting the projected density of states, this object is also used to exctract the band edges at the left electrode. These informations are printed and used to align the electrostatic difference potential plots.
  • ElectrostaticDifferencePotential: to plot the average electrostatic difference potential across the device.

In order to get the DeviceDensityOfStates plot without an applied bias, run the script from a terminal:


where contains the analysis objects described above. You can also create plots for finite bias. See the section Finite-bias calculations. Results are illustrated below.

../../_images/DDOS_0bias.png ../../_images/DDOS_minus1bias.png ../../_images/DDOS_1bias.png

Finite-bias calculations

In order to create the finite-bias plots (Vbias = -1 and 1 V), do as follows:

  1. Go to the Labfloor and from the file select the device configuration obtained at Vbias = -1 V, and drag and drop it onto the script_generator_icon Scripter. The following window will show up.

  2. Remove the calculator_icon New Calculator object and add instead analysis_from_file_icon Analysis from File.

  3. In analysis_from_file_icon Analysis from File, set the NetCDF file and the Object id of the device configuration object.

  4. Then, add the Electrondensity, DeviceDensityOfStates and ElectrostaticDifferencePotential analysis objects and set the parameters as described above.

  5. Change the output file to

  6. Before running the job, send the script to the editor_icon Editor, and add the following line to get the device_configuration file within the output file.

    # -------------------------------------------------------------
    # Analysis from File
    # -------------------------------------------------------------
    device_configuration = nlread('', object_id='ivcurve010')[0]
  7. After this part and just before ElectronDensity, add the following line:

    nlsave('', device_configuration)
  8. Lastly, save the script and run the job. Once the job have finished, run the in order to get the Device Density of States plot. Run the script as follow:


You can plot the Device Density of States plot following the same procedure for the ivcurve_selftconsistent_configuration.ncivcurve020 file.

Electrostatic potential


  1. In the Labfloor, select the ElectrostaticDifferencePotential objects obtained for Vbias = 0, -1, 1,

  2. Use the 1D Projector plugin to plot the data.


You can see that the results are not perfectly converged wrt. the length of the device, because the gradient of the electrostatic potential is not zero at the ends of the device central region. The figure below compares the ElectrostaticDifferencePotential of the device used in this tutorial with results for a longer one.



[aTB09]F. Tran and P. Blaha. Accurate band gaps of semiconductors and insulators with a semilocal exchange-correlation potential. Phys. Rev. Lett., 102:226401, 2009. doi:10.1103/PhysRevLett.102.226401.