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¶
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.
You will now set up an ATK semiempirical (ATK-SE) calculation using the Slater-Koster model Hamiltonian.
Add the following blocks to the Script by double-clicking the correspoding icons in the Blocks panel:
- 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.
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.
You can zoom into a region close to the band edges and precisely measure the calculated indirect band gap, 1.186 eV.
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.
- 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.
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 using the 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):
On the LabFloor you will now have four different Bandstructure objects included in the “MGGA_bulk.nc”.
You will use the following script to fit the c-parameter: bandstructure_fit.py. Save it in the project directory.
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 “bandstructure_fit.py”, you will see that band gaps are calculated from a Bandstructure object like this:
bandstructure = nlread("file.nc", Bandstructure) 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
You will now set up the silicon device.
- Navigate to .
- In the Surface (Cleave) window, choose the [1,0,0] Miller indices and click “Next”.
Keep the 1\(\times\)1 surface unit cell and click “Next”.
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.
Your Si(100) configuration is now constructed, and you will use it to construct a device with the transport direction along (100). Navigate to, keep the default parameters, and simply click “OK”.
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 .
- Type in c < 0.5 and press Enter to select all atoms in the left hand side of the device.
- Use 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.
- 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.
You can now add all the analysis objects needed. Add the following ones:
- Increase considerably the k-point sampling. Use a 21x21 grid.
- Energy: from -2 to 2 eV with 401 points.
Change the output filename to something like
Doping the silicon junction
device_configuration = DeviceConfiguration( central_region, [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) ]) central_region.setExternalPotential(external_potential) # Add doping to left and right electrodes left_electrode.setExternalPotential( AtomicCompensationCharge([(Silicon, -1*doping)])) right_electrode.setExternalPotential( AtomicCompensationCharge([(Silicon, doping)]))
Running 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.
DFT-MGGA IV curve¶
- Choose ATK-DFT.
- Exchange correlation: MGGA.
- SingleZetaPolarized basis set.
Analyzing the results¶
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:
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.
In order to compare the Slater-Koster and MGGA IV curves, download the script IV_compare.py to your project folder and send it to the 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 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.
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: ddos_edp.py.
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 and customize the plot according to your wish.
In particular, the ddos_edp.py 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:
atkpython ddos_edp.py IV_2e19_SK.nc
IV_2e19_SK.nc contains the analysis objects described above. You can also create plots for finite bias. See the section Finite-bias calculations. Results are illustrated below.
In order to create the finite-bias plots (Vbias = -1 and 1 V), do as follows:
Go to the Labfloor and from the
ivcurve_selftconsistent_configuration.ncfile select the device configuration obtained at Vbias = -1 V, and drag and drop it onto the Scripter. The following window will show up.
Then, add the Electrondensity, DeviceDensityOfStates and ElectrostaticDifferencePotential analysis objects and set the parameters as described above.
Change the output file to
# ------------------------------------------------------------- # Analysis from File # ------------------------------------------------------------- device_configuration = nlread('ivcurve_selfconsistent_configurations.nc', object_id='ivcurve010')
After this part and just before ElectronDensity, add the following line:
Lastly, save the script and run the job. Once the job have finished, run the ddos_edp.py in order to get the Device Density of States plot. Run the ddos_edp.py script as follow:
atkpython ddos_edp.py Analysis_ivcurve_10.nc
You can plot the Device Density of States plot following the same procedure for the
In the Labfloor, select the ElectrostaticDifferencePotential objects obtained for Vbias = 0, -1, 1,
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.