# 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](../figs/si_soi_tilted.png)
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](../figs/si_soi_regions.png)

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.

```python
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 = cal.calc_charge_current()
print(I)

```

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

```python
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,

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

```python
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.

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

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

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

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

```python
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](../figs/si_soi_work_func.png)

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

```python
#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.

```python
from nanotools import TwoProbe
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](../figs/si_soi_vdh.png)

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](../figs/si_soi_iv.png)

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

```python
from nanotools import TwoProbe, 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`,

```python
from nanotools 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

```python
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](../figs/si_soi_field_avg.png)

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