Silicon nanowire field-effect transistor¶
This tutorial shows you how to set up and perform calculations for a device based on a silicon nanowire. You will define the structure of a hydrogen passivated Si(100) nanowire, and set up a field-effect transistor (FET) structure with a cylindrical wrap-around gate. The tutorial ends with calculations of the device conductance as a function of gate bias.
You will primarily use the graphical user interface QuantumATK for setting up and analyzing the results. If you are not familiar with QuantumATK, please go through the tutorial Introduction to QuantumATK.
The underlying calculation engines for this tutorial are ATK-DFT and ATK-SE. A complete description of all the parameters, and in many cases a longer discussion about their physical relevance, can be found in the ATK Reference Manual.
Band structure of a Si(100) nanowire¶
First step is to set up the Si(100) nanowire and optimize the geometry. You should use the ATK-DFT calculator for this. We then compute the band structure of the nanowire using 3 different computational models; DFT-GGA, DFT-MetaGGA, and the Extended Hückel method.
Building the nanowire¶
- Go to and type “silicon fcc” to locate the diamond phase of silicon.
- Add the Silicon (alpha) bulk configuration to the Stash (double-click or use the icon).
In the Builder, use thetool to create the Si(100) facet:
- keep the default (100) cleave direction, and click Next;
- keep the default surface lattice, and click Next;
- keep the default supercell, which will ensure that the wire direction is perpendicular to the surface, and click Next.
- Click Fińish to add the cleaved structure to the Stash.
Next, use thetool to repeat the structure twice along the A and B directions:
Finalize the nanowire by following these steps:
- use the tool to set the A and B lattice vector lengths to 20 Å;
- use the tool to center the structure along all directions;
- go to and apply the 4(SP3) type hydrogen passivation to the nanowire.
Setting up and running the calculations¶
You should now relax the nanowire and compute the band structure. The DFT-GGA method is used for geometry optimization, while the band structure is computed using DFT-GGA, DFT-MetaGGA, and the Extended Hückel model.
The Meta-GGA and Extended Hückel models can not be used for relaxation.
- New Calculator
- New Calculator
- New Calculator
- set the k-point sampling to 1x1x11;
- change the exchange-correlation potential to GGA.
- set the k-point sampling to 1x1x11;
- change the exchange-correlation potential to MGGA.
- select the ATK-SE: Extended Hückel calculator;
- uncheck “No SCF iteration” to make the calculation selfconsistent;
- set the k-point sampling to 1x1x11;
- increase the density mesh cut-off to 20 Hartree;
- go to the Hückel basis set tab, and select the “Hoffmann.Hydrogen” and “Cerda.Silicon [GW diamond]” basis sets. The latter has been fitted to GW calculations, and gives an excellent description of the silicon bandstructure, including the value of the band gap.
The default value for the density mesh cut-off is often sufficient, but by increasing it you may sometimes obtain a slightly higher accuracy at a small additional computational cost.
You should also adjust the value of the c-parameter used in the TB09 MetaGGA method. This is essentially a fitting parameter that can be used to adjust the calculated band gap of the material. If the c-parameter is not defined by the user, a value is automatically calculated from the electronic structure of the material. However, this will not work for a system with vacuum regions, like a nanowire.
You will therefore use a fixed value, c=1.0364, which has been fitted to yield an accurate prediction of the fundamental band gap in silicon (1.13 eV).
Transfer the QuantumATK Python script to the Editor using the icon.
Then locate the script line defining the variable
exchange_correlation for the MGGA
exchange-correlation, and set the fixed value of c to 1.0364, as illustrated below.
#---------------------------------------- # Exchange-Correlation #---------------------------------------- exchange_correlation = MGGA.TB09LDA(c=1.0364)
Running the calculations¶
Analyzing the results¶
Once the calculations have finished, the NetCDF data file
should appear on the QuantumATK LabFloor. Locate the Bandstructure items inside this file.
The order of the item IDs corresponds to the order in which the band structures were
computed and saved, i.e. GGA, MGGA, and Hückel.
Select one of the band structures and open the Bandstructure Analyzer in the right-hand Panel Bar to inspect the computed band structure. You can do this for all three band structures, or use the Compare Data plugin to plot a direct comparison of the three band structures.
Detailed band structure analysis¶
In order to compare the band gaps obtained with the three computational models,
it is convenient to use
bandstructure_analyzer.py, which uses QuantumATK Python for detailed analysis
of the gaps in a band structure. The script is shown below:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141
from QuantumATK import * import pylab # Custom analyzer for calculating the band gap of a bandstructure #helper function to find minima def fitEnergyMinimum (i_min, energies, k_points, nfit=3): """ Function for fitting the energy minimum located around i_min. @param i_min : approximate position of the energy minimum. @param energies: list of energies. @param k_points: list of k_points which correspond to the energies. @param nfit : order of the polynomium. @return d2_e, e_min, k_min : second derivative of energy, minimum energy, minimum k_point. """ #list of energies n = len(energies) efit = numpy.array([energies[(n+i_min-nfit/2+i)%n] for i in range(nfit)]) kfit = numpy.array([i-nfit/2 for i in range(nfit)]) #special cases if i_min == 0: #assume bandstructure symmetric around zero for i in range(nfit/2): efit[i] = energies[nfit/2-i] if i_min == n-1: #assume bandstructure symmetric around end point for i in range(nfit/2+1, nfit): efit[i] = energies[n-1+nfit/2-i] #make fit p = numpy.polyfit(kfit,efit,2) i_fit_min = -p/2./p pf = numpy.poly1d(p) e_min = pf(i_fit_min) i0 = int(i_fit_min+i_min+n) w = i_fit_min+i_min+n-i0 k_min = (1-w)*k_points[i0%n]+w*k_points[(i0+1)%n] return p, e_min, k_min def analyseBandstructure(bandstructure, spin): """ Function for analysing a band structure and calculating bandgaps. @param bandstructure : The bandstructure to analyze @param spin : Which spin to select from the bandstructure. @return e_val, e_con, e_gap : maximum valence band energy, minimum conduction band energy, and direct band gap. """ energies = bandstructure.evaluate(spin=spin).inUnitsOf(eV) #some dummy variable to help finding the extrema e_valence_max = -1.e10 e_conduction_min = 1.e10 e_gap_min = 1.e10 i_valence_max = 0 i_conduction_min = 0 i_gap_min = 0 n_valence_max = 0 n_conduction_min = 0 n_gap_min = 0 # Locate extrema for i in range(energies.shape): # find first state below Fermi level n = 0 while n < energies.shape and energies[i][n] < 0.0: n += 1 # find maximum of valence band if (energies[i][n-1] > e_valence_max): e_valence_max = energies[i][n-1] i_valence_max = i n_valence_max = n-1 # find minimum of conduction band if (energies[i][n] < e_conduction_min): e_conduction_min=energies[i][n] i_conduction_min=i n_conduction_min=n # find minimum band gap if (energies[i][n]-energies[i][n-1] < e_gap_min): e_gap_min = energies[i][n]-energies[i][n-1] i_gap_min = i n_gap_min = n-1 # Print out results a_val, e_val, k_val = fitEnergyMinimum(i_valence_max, energies[:,n_valence_max], bandstructure.kpoints()) print 'Valence band maximum %7.4f eV at [%6.4f, %6.4f,%6.4f] ' \ %(e_val, k_val, k_val, k_val) a_con, e_con, k_con = fitEnergyMinimum(i_conduction_min, energies[:,n_conduction_min], bandstructure.kpoints()) print 'Conduction band minimum %7.4f eV at [%6.4f, %6.4f,%6.4f] ' \ %(e_con, k_con, k_con, k_con) print 'Fundamental band gap %7.4f eV ' % (e_con-e_val) a_gap, e_gap, k_gap = fitEnergyMinimum(i_gap_min, energies[:,n_gap_min+1]- energies[:,n_gap_min], bandstructure.kpoints()) print 'Direct band gap %7.4f eV at [%6.4f, %6.4f,%6.4f] ' \ %(e_gap, k_gap, k_gap, k_gap) return e_val, e_con, e_gap def analyzer(filename, **args): """ Find band gaps of band structures in netcdf file. """ if filename == None: return #read in the bandstructure you would like to analyze bandstructure_list = nlread(filename, Bandstructure) if len(bandstructure_list) == 0 : print 'No Bandstructures in file ', filename return for s in [Spin.All]: b_list =  n = 0 for b in bandstructure_list: print 'Analyzing bandstructure number ', n b_list = b_list + [analyseBandstructure(b,s)] print print n += 1 x = numpy.arange(len(b_list)) e_val = numpy.array([b for b in b_list]) e_con = numpy.array([b for b in b_list]) e_indirect = e_con-e_val e_direct = numpy.array([b for b in b_list]) analyzer("si_100_nanowire.nc")
Download the script (link:
bandstructure_analyzer.py), and execute it using the Job Manager.
You should get the following output:
Analyzing bandstructure number 0 Valence band maximum -1.6395 eV at [0.0000, 0.0000,0.0121] Conduction band minimum 1.6400 eV at [0.0000, 0.0000,0.0000] Fundamental band gap 3.2795 eV Direct band gap 3.2799 eV at [0.0000, 0.0000,0.0000] Analyzing bandstructure number 1 Valence band maximum -1.8434 eV at [0.0000, 0.0000,0.1501] Conduction band minimum 1.8450 eV at [0.0000, 0.0000,0.0000] Fundamental band gap 3.6884 eV Direct band gap 3.6914 eV at [0.0000, 0.0000,0.0026] Analyzing bandstructure number 2 Valence band maximum -1.7737 eV at [0.0000, 0.0000,0.0000] Conduction band minimum 1.7735 eV at [0.0000, 0.0000,0.0000] Fundamental band gap 3.5472 eV Direct band gap 3.5472 eV at [0.0000, 0.0000,0.0000]
The output shown above gives the band structure analysis for the DFT-GGA, DFT-MGGA, and Extended Hückel models. The script gives the valence band maximum, conduction band minimum, fundamental band gap, and direct band gap from the valence band maximum. The fundamental band gap for the Si(100) nanowire is significantly larger than that of bulk silicon, and all three computational models place it at the Gamma point or close to it.
Furthermore, the MetaGGA model gives a gap which is about 0.5 eV larger than the GGA model, which is in accordance with the difference in the band gaps of bulk Si calculated with the GGA and MetaGGA methods. The Extended Hückel model is in agreement with the MetaGGA method.
Si(100) nanowire FET device¶
You will now consider a Si(100) nanowire field-effect transistor with a cylindrical wrap-around gate. You should first build the device configuration, including the gate and doping of the silicon nanowire, and then use the Extended Hückel model to compute device properties such as zero-bias transmission and conductance vs. gate voltage.
FET device configuration¶
You should use the relaxed nanowire configuration as a starting point for building
the device. The contents of the file
si_100_nanowire.nc should be available
on the QuantumATK LabFloor. It contains four bulk configurations: The first one (gID000)
is the unrelaxed structure, while the remaining ones correspond to the relaxed structure,
calculated with different methods (DFT-GGA, DFT-MGGA and Extended Hückel).
Setting up the gate¶
Next, use thetool to define the metallic wrap-around gate:
- First, add a new metallic region with a value of 0 Volt.
- Under Geometry, select Tube to create a cylindrical region.
- Define the geometry of the tube by entering the parameters shown in the image below. The cylinder will extend to the edge of the simulation cell along A and B, and will cover most of the central part of the nanowire along C.
Build the device¶
Use thetool to build the Si(100) nanowire FET device from the bulk nanowire configuration. Use the default suggestion for the electrode lengths.
Doping the Si(100) wire¶
You will now introduce doping in the nanowire. Instead of explicitly adding dopant atoms, which would result in a very high doping concentration for this relatively small device, you will set up the system as a p-i-n junction, by adding certain amounts of charge to the electrodes. This will lead to shifts in the Fermi levels of the two electrodes, resulting in a built-in potential in the device.
Use the mouse to select all atoms in the left electrode and open thetool. In the window that shows up, choose the following settings:
- Doping Type: p-type
- Value: 4e+19
- Unit: e/cm3 (left electrode)
Next, close the Doping window and select all atoms in the right electrode. Open again the Doping tool, and choose the following settings:
- Doping Type: n-type
- Value: 4e+19
- Unit: e/cm3 (right electrode)
Zero gate voltage calculation¶
You will now calculate the transmission spectrum of the p-i-n doped Si(100) nanowire
FET device at zero gate potential using the Extended Hückel model. In the
Script Generator, change the default output file name to
and add the following blocks to the script:
- Uncheck “No SCF iteration” to make the calculation selfconsistent.
- Increase the density mesh cut-off to 20 Hartree.
- Go to the Hückel basis set tab and select the “Hoffmann.Hydrogen” and “Cerda.Silicon [GW diamond]” basis sets.
- Under Poisson solver, use Neumann boundary conditions in the A and B directions.
In the TransmissionSpectrum block, set the energy range to -4 to +4 eV with 301 points.
It would be physically incorrect to use periodic boundary conditions along the A and B directions when there are gates present in the system. We therefore choose Neumann boundary conditions along A and B instead.
Analyzing the results¶
Once the calculation is completed, locate the
NetCDF file on the Labfloor and unfold it. Select the TransmissionSpectrum item
and use the Transmission Analyzer to plot the transmission spectrum.
The transmission spectrum is shown above. The transmission is zero inside the nanowire band gap, which appears to be roughly 7 eV, significantly larger than computed earlier. This is due to the p- and n-doping, which has moved the valence and conduction band edges with respect to the Fermi level.
Performing a gate scan¶
Finally, you will calculate the transmission spectrum for different values of the gate bias. You can then plot the device conductance as function of gate bias.
The calculations are most conveniently done using QuantumATK Python scripting. Use the script
nanodevice_gatescan.py for this. Download it to the active QuantumATK project folder, which also contains
si_100_nanowire_fet_pin.nc. The script is reproduced below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39
# Read in the old configuration device_configuration = nlread("si_100_nanowire_fet_pin.nc",DeviceConfiguration) calculator = device_configuration.calculator() metallic_regions = device_configuration.metallicRegions() # Define gate_voltages gate_voltage_list=[-1, 0, 1.0, 2.0, 3.0, 4.0, 5.0, 6, 7, 8,]*Volt # Define output file name filename= "nanodevice_huckel.nc" # Perform loop over gate voltages for gate_voltage in gate_voltage_list: # Set the gate voltages to the new values new_regions = [m(value = gate_voltage) for m in metallic_regions] device_configuration.setMetallicRegions(new_regions) # Make a copy of the calculator and attach it to the configuration # Restart from the previous scf state device_configuration.setCalculator(calculator(), initial_state=device_configuration) device_configuration.update() nlsave(filename, device_configuration) # Calculate analysis objects electron_density = ElectronDifferenceDensity(device_configuration) nlsave(filename, electron_density) electrostatic_potential = ElectrostaticDifferencePotential(device_configuration) nlsave(filename, electrostatic_potential) transmission_spectrum = TransmissionSpectrum( configuration=device_configuration, energies=numpy.linspace(-4,4,301)*eV, kpoints=MonkhorstPackGrid(1,1), energy_zero_parameter=AverageFermiLevel, infinitesimal=1e-06*eV, ) nlsave(filename, transmission_spectrum) nlprint(transmission_spectrum)
Execute the script using the Job Manager or from a command
line. The script consists of 8 different calculations and it will take about 5 hours
running in parallel with 4 MPI processes. You can, however, decide to reduce this
time by reducing the number of gate voltages that the transmission spectrum is computed
for – simply edit the variable
Analyzing the gate scan¶
Download and execute the script
conductance_plot.py, which reads the computed transmission spectra,
calculates the conductance for each, and plots the conductance against gate bias.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28
# Read the data transmission_spectrum_list = nlread("nanodevice_huckel.nc", TransmissionSpectrum) configuration_list = nlread("nanodevice_huckel.nc", DeviceConfiguration) conductance = numpy.zeros(len(configuration_list)) gate_bias = numpy.zeros(len(configuration_list)) for i, configuration in enumerate(configuration_list): transmission_spectrum = transmission_spectrum_list[i] energies = transmission_spectrum.energies().inUnitsOf(eV) spectrum = transmission_spectrum.evaluate() gate_bias[i] = configuration.metallicRegions().value().inUnitsOf(Volt) conductance[i] = transmission_spectrum.conductance() # Sort the data according to the gate bias index_list = numpy.argsort(gate_bias) # Plot the spectra import pylab pylab.figure() ax = pylab.subplot(111) ax.semilogy(gate_bias[index_list], conductance[index_list]) ax.set_ylabel("Conductance (S)", size=16) ax.set_xlabel("Gate Bias (Volt)", size=16) ax.set_ybound(lower=1e-26) for g,c in zip(gate_bias[index_list], conductance[index_list]): print g,c pylab.show()