Silicon Transistor tutorial

Requirements

Software components

  • nanotools

  • NanoDCAL+

Structure file

  • slab.xyz

Pseudopotentials

  • H_AtomicData.mat

  • Si_AtomicData.mat

Copy the mat files from the pseudo database to the current working directory and export NANODCALPLUS_PSEUDO=$PWD.

References

Briefing

In this tutorial, we demonstrate simulating the on-off transport property of a npn silicon based junction. We apply a source-drain bias over the junction and compute the charge current for a series of gate voltages. The atomic structure of the transport junction is shown in the following figure transmission profile of carbon nano-tube The z-direction is chosen as the transport direction.

The whole system has a slab geometry, periodic in the x-direction and passivated by hydrogen atoms in the y-direction. Above and below the passivated surfaces are dielectric media as indicated in the following figure

transmission profile of carbon nano-tube

P-doped silicon atoms are highlighted in green, and the rest are n-doped. We introduce two gates at the top and bottom of the simulation box, so as to tune the electrostatic potential within the transport channel; without gates, the junction is naturally in the off-state.

Full script

The following script is used to compute the charge current under a gate-source voltage of 0.5 V and a drain-source bias of 0.2 V. The result, a single number with proper unit, is printed out on the screen at the end. Besides, standard output files ending with _out.json or _out.h5 are saved in the directory after the program is finished.

import numpy as np
from nanotools import Atoms, Cell, Current, Pop, Region, System, TotalEnergy, TwoProbe, IV
npc = 40    # number of processors to be used

cell = Cell(avec=[[5.43, 0.0, 0.0], [0.0, 30., 0.0], [0.0 , 0.0 , 5.43]], resolution=0.2)
atoms= Atoms(positions="slab.xyz", doping_ratios="slab.xyz")

pop = Pop(type="fd", contour_type="single", sigma=0.026)
sys = System(cell=cell, atoms=atoms, pop=pop)
sys.kpoint.set_grid([5,1,21])

di1 = Region(fractional_x_range = [0,1], fractional_y_range = [0,0.25], fractional_z_range = [0,1])
di2 = Region(fractional_x_range = [0,1], fractional_y_range = [0.75,1], fractional_z_range = [0,1])
sys.add_dielectric(di1, epsilon=25)
sys.add_dielectric(di2, epsilon=25)
sys.set_open_boundary(['-y','+y'])
left = TotalEnergy(sys)
left.solver.set_mpi_command(f"mpiexec -n 4")
left.solver.restart.densityPath = "left_out.h5"

center = left.copy()
center.supercell(T=np.diag([1,1,30]))
center.system.kpoint.set_grid([5,1,1])
dopp = [.999, .999, .999, 1, 1, .999, .999, .999, .999, 1, 1,
        .999, .999, .999, 1, 1, .999, .999, .999, .999, 1, 1]
dopn = [1.001, 1.001, 1.001, 1, 1, 1.001, 1.001, 1.001, 1.001, 1, 1,
        1.001, 1.001, 1.001, 1, 1, 1.001, 1.001, 1.001, 1.001, 1, 1]
center.system.atoms.set_doping_ratios(dopn*10 + dopp*10 + dopn*10)

g1  = Region(fractional_x_range = [0,1], fractional_y_range = [0,0], fractional_z_range = [.33,.66])
g2  = Region(fractional_x_range = [0,1], fractional_y_range = [1,1], fractional_z_range = [.33,.66])
vgs = +0.5  # gate-source voltage Vg-Vs
wf  = 4.5   # work function
center.system.add_gate(region=g1, work_func=wf, Vgs=vgs)
center.system.add_gate(region=g2, work_func=wf, Vgs=vgs)

center.solver.mix.alpha = 0.1
#center.solver.mix.precond = "kerker"
center.solver.mix.maxit = 120
center.solver.set_mpi_command(f"mpiexec -n {npc}")
center.solver.mix.tol = 1.0e-6
center.solver.mpidist.kptprc = 5  # k-parallelization setting
center.solver.mpidist.zptprc = 4
center.solver.restart.densityPath = "center_out.h5"

vds = 0.2   # drain-source voltage Vd-Vs
dev = TwoProbe(left=left, center=center, transport_axis=2, bias=vds)
dev.solve() # self-consistent calculation

dev.center.system.kpoint.set_grid([25,1,1])  # increase k-points for current computation
cal = Current.from_twoprobe(dev)
I,U = cal.calc_charge_current()
print(I)
print(U)

Note that it may take over 60 iterations to reach convergence. In that case, the script needs be run several times depending on center.solver.mix.maxit parameter. Carefully tune center.solver.mix.alpha if it helps.

Explanations

Doping level

Just as most semiconductor devices, our silicon-based transistor requires a certain doping level in order to be functionalized. In our code, the doping feature is implemented in an artificial manner (virtual crystal approximation), namely no actual dopants are added into the system. Instead, we specify a doping ratio on each individual atom in the system. A ratio of 1.001 (0.999) means 0.1% of the atom’s valance electrons are to be added (removed).

When we initialize the atoms object as follows

atoms= Atoms(positions="slab.xyz", doping_ratios="slab.xyz")

we set the doping_ratios keyword to be the name of the atomic structure file “slab.xyz”. This means that the doping ratio on each atom is specified in the structure file. Taking a look into the file,

22
AtomType   X               Y            Z         doping_ratio
Si        2.98168995  12.301750000  0.26668995    1.001
Si        0.26668995  15.016750000  0.26668995    1.001
Si        2.98168995  17.731750000  0.26668995    1.001
H         0.80006984  10.000000000  0.80006984    1
H         5.16331005  20.033523324  0.80006984    1
...

we see that the ratios are listed under the doping_ratio keyword.

It is also possible to set doping ratios for an existing Atoms object, using its set_doping_ratios method. This method takes in a list of numbers with the same order as the atomic coordinates. For example, to set the doping ratios in the central region, which consists of ten n-doped unit cells, followed by 10 p-doped cells, and then followed by another ten n-doped unit cells, we use

dopn = [1.001, 1.001, 1.001, 1, 1, 1.001, 1.001, 1.001, 1.001, 1, 1,
        1.001, 1.001, 1.001, 1, 1, 1.001, 1.001, 1.001, 1.001, 1, 1]
dopp = [.999, .999, .999, 1, 1, .999, .999, .999, .999, 1, 1,
        .999, .999, .999, 1, 1, .999, .999, .999, .999, 1, 1]
center.system.atoms.set_doping_ratios(dopn*10 + dopp*10 + dopn*10)

Remember that the “*” operator for lists means replication and that “+” means concatenation.

Dielectrics and gates

In order to define dielectrics and gates within the system, we first need to create some Region objects which specify the geometry of the desired object.

di1 = Region(fractional_x_range = [0,1], fractional_y_range = [0,0.25], fractional_z_range = [0,1])
di2 = Region(fractional_x_range = [0,1], fractional_y_range = [0.75,1], fractional_z_range = [0,1])

Region object specifies a rectangular/parallelepiped region in the 3d simulation box. The range in each dimension is specified by a pair of fractional numbers. Once a system object is initialized, a dielectric region can be added using the add_dielectric method with the specified (relative) permittivity, denoted as epsilon, eg

system.add_dielectric(di1, epsilon=25)
system.add_dielectric(di2, epsilon=25)

It is important to specify the proper boundary condition of Poisson equation in accordance with the dielectric configuration. In our case, we put

system.set_open_boundary(['-y','+y'])

which indicates that the Neumann boundary condition is applied to the surfaces normal to y-axis.

Similarly we may create regions for metallic gates

g1  = Region(fractional_x_range = [0,1], fractional_y_range = [0,0], fractional_z_range = [.33,.66])
g2  = Region(fractional_x_range = [0,1], fractional_y_range = [1,1], fractional_z_range = [.33,.66])

and add gates to the system object

system.add_gate(region=g1, work_func=wf, Vgs=vgs)
system.add_gate(region=g2, work_func=wf, Vgs=vgs)

where Vgs denotes the gate-source voltage, and work_func denotes the work-function of the gate material. Note that, unlike dielectric regions, for gates only 2d geometry is supported. That is why we collapse the fractional_y_range to [0,0] or [1,1], i.e. only on the surfaces of the simulation box. The effect of gate is providing a Dirichlet boundary condition for Poisson equation. Only one gate is allowed on each surface, and Neumann boundary condition is always assumed on the rest of the surface apart from the gate.

The gate work-function is a parameter we expect to be provided by the user, especially when the user has a specific gate material in mind. Otherwise, we refer to the next section as for estimating a proper work function in case the user simply wants a virtual gate.

Work function of gates

Adding metallic gates to the surface of our simulation box essentially sets a boundary value for the electrostatic potential in the system. However, it is important to take note that this boundary value differs from the gate voltage we apply. In fact, it should be the gate’s work function minus the gate-source voltage. See the diagram below.

transmission profile of carbon nano-tube

The diagram depicts an imaginary contact between the source and gate materials. As a convention, the chemical potential (Fermi energy, red dashed line) of source is set at zero; that of the gate is somewhat lower as we apply a positive gate-source voltage (remembering electrons have a negative charge). The boundary value provided by the gate is actually the vacuum level of the material, which is equal to its work function minus the gate-source voltage.

Having understood the role of work function, we now demonstrate how to calibrate it in case the user does not have this material parameter at hand or simply wants an effective way to tune the electrostatics via a virtual gate. The goal of this calibration is such that the boundary value stays approximatly the same, whether without the gate at all or with the gate in equilbrium to the source (Vgs=0). To this end, we comment out

#center.system.add_gate(region=g1, work_func=wf, Vgs=vgs)
#center.system.add_gate(region=g2, work_func=wf, Vgs=vgs)

and set the drain-source bias Vds to zero. After the self-consistent calculation is finished, run the following script to visualize the electrostatic potential on the surface where the gate is to be placed.

ecalc = TwoProbe.read(filename="nano_2prb_out.json")
veff = ecalc.smooth_field("potential/effective", axis=2, preslice=lambda a: a[:,-2:-1,:])
fig = ecalc.center.plot_field(veff, axis=2)

transmission profile of carbon nano-tube

One can see that the electrostatic potential is bumped up in the middle due to the p-n junction effect. The potential in the middle reads 4.5 V. That is what we assign to the gate work function.

Results

Run the script over and over under a series of gate-source voltages, and read off the current values printed on screen at the end. The resulting I-Vgs curve looks like transmission profile of carbon nano-tube

We have implemented an IV calculator to automate this process, i.e. run self-consistent calculations under a series of voltages, compute the charge currents, collect data and plot the I-V curve. The usage is as follows

from nanotools.iv import IV
# see the script at the top ...
dev = TwoProbe(left=left, center=center, transport_axis=2, bias=vds)
calc = IV(reference_calculator=dev,
          k_grid_4_current=[25,1,1],
          varying_elements=center.system.hamiltonian.extpot.gates,
          voltages=[0.0,0.3,0.6,0.9,1.0,1.1,1.2],
          work_func=wf)
calc.solve()
calc.plot(log=True, filename="iv.png")

Or, if there is already an existing nano_2prb_in.json,

from nanotools.iv import TwoProbe, IV
dev = TwoProbe.read("nano_2prb_in.json")
calc = IV(reference_calculator=dev,
          k_grid_4_current=[25,1,1],
          varying_elements=dev.center.system.hamiltonian.extpot.gates,
          voltages=[0.0,0.3,0.6,0.9,1.0,1.1,1.2],
          work_func=4.5)
calc.solve()
calc.plot(log=True, filename="iv.png")

The IV calculator is initialized with the reference TwoProbe calculator, a varying_elements list containing the gate objects whose voltages we wish to vary (in our case we wish to vary both gates simultaneously), a range of voltages in voltages, and the work-function of the gates. calc.solve() then carries out self-consistent calculations and computes the corresponding currents, which are stored in calc.currents attribute. Finally, calc.plot(log=True) plots the I-V curve in log-scale.

Note that, during the process, there is no guarantee that the self-consistent calculations are well converged: they may stop simply because the number of iterations has hit mixer.maxit. In this case, the user may want to run calc.solve() several times until convergence is reached under each voltage. Alternatively, the user can enter each of the spawned directories, whose name ends with the corresponding voltage value, and manually launch computations from therein.

From the above figure, we see that the current gets switched on as the gate voltage goes more positive, as expected for a npn junction. To have a closer look, we plot the effective potential using the following script

ecalc = TwoProbe.read(filename="nano_2prb_out.json")
veff = ecalc.smooth_field("potential/effective", axis=2, preslice=lambda a: a[:,60:90,:])
fig = ecalc.center.plot_field(veff, axis=2, filename="potential.png")

on the results obtained under Vgs=0 V and 1 V respectively. The figure is shown below

transmission profile of carbon nano-tube

As one can see, the potential barrier for electrons gets significantly lowered at Vgs=1 V, hence allowing for more charge flow.