# 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 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

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.

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)
```

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

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

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