5. Poisson–NEGF–Master-equation simulations of an FD-SOI field-effect transistor with a single quantum dot

5.1. Requirements

5.1.1. Software components

  • QTCAD

  • Gmsh

5.1.2. Geometry file

  • qtcad/examples/tutorials/meshes/fdsoi_poisson_negf_master.geo

5.1.3. Python script

  • qtcad/examples/tutorials/poisson_negf_master.py

5.1.4. References

5.2. Briefing

In this tutorial, we will compute the Coulomb peaks and the charge stability diagram of a single quantum dot in a fully-depleted silicon-on-insulator (FD-SOI) structure using a computationally efficient, yet accurate workflow. In this workflow, the device electrostatics are resolved using a simple non-linear Poisson equation. The current through the device in the single-electron, near-equilibrium regime is then calculated using a computationally-light calculation based on the nonequilibrium Green’s function (NEGF) formalism, which overcomes several limitations of the WKB equation introduced in Quantum transport—WKB approximation; see, for example, our previous tutorial: Quantum transport—Master equation. The result of the single-electron NEGF calculation is then used to calibrate the current predicted by the master-equation solver, which is itself then used to calculate the current in the many-electron regime. Finally, the Coulomb peaks and the charge stability diagram are obtained.

Overall, this workflow combines the accuracy of single-electron NEGF calculations with the computational efficiency of the master equation solver to obtain an accurate description of charge transport in the many-electron regime.

5.3. Computational approach

In the first step, we solve the electrostatics of the device using the non-linear Poisson equation, which assumes thermodynamic equilibrium. For weak source–drain bias, the device is sufficiently close to equilibrium to use the electrostatics obtained from the Poisson equation to calculate the current through the channel and the leads within the NEGF formalism.

Because the NEGF approach models quantum transport using a single-body Hamiltonian, it does not accurately describe current in the Coulomb blockade regime, when multiple interacting electrons are present. In contrast, the master equation solver uses the many-body Hamiltonian to accurately model Coulomb blockade. Therefore, we use the NEGF formalism to calculate the current for single electrons in the sequential tunneling regime. More specifically, using the NEGF formalism, we calculate the current associated with the first Coulomb peak, which corresponds to tunneling of a single electron through an empty quantum dot.

In the final step, we solve the master equation to obtain the current in the few-electron regime using the featureless approximation, while globally rescaling all current values in such a way that the amplitude of the first peak matches the current obtained in the NEGF simulation.

The approximations made to compute the current through the device in this tutorial are described below.

  • Lever-arm approximation: As described in Generalization to multiple-dot systems: Lever arm matrix, we may approximate the effect of plunger gate bias variations as linear changes in the quantum-dot single-particle energy levels through a lever-arm matrix, which is closely related to the capacitance matrix of the system. Through this approximation, one can avoid solving the Poisson, single-particle Schrödinger, and many-body Schrödinger equations again for each new bias configuration when computing charge stability diagrams.

  • Near-equilibrium approximation: In this tutorial, we calculate the electrostatics at equilibrium and use the corresponding conduction band edge to compute the source–drain current in the NEGF formalism.

  • Constant broadening function: Finally, we use a constant broadening function in the master-equation solver to obtain the current in the featureless approximation. This assumes that the broadening function is approximately independent of energy over the source–drain bias range, in addition to being approximately independent of the single-particle basis state index.

5.4. Device geometry and mesh construction

The device considered in this tutorial is very similar to the single-dot FD-SOI FET presented in NEGF–Poisson simulations of an FD-SOI field-effect transistor with a quantum dot. It consists of a silicon nanowire with a rectangular cross section which is surrounded by silicon dioxide. The silicon channel has a width of \(10~\textrm{nm}\) and a thickness of \(6~\textrm{nm}\). It is partitioned into four regions: a degenerately-n-doped source, an intrinsic channel, a quantum dot region, and a degenerately-n-doped drain.

To take advantage of our implementation of the NEGF formalism, the channel and lead regions of the device must be meshed in a structured manner, in which the mesh contains well-defined principal layers along the transport direction. This requirement can be achieved by either using the extrusion or transfinite meshing methods available in Gmsh. In this tutorial, we start from the geometry file and then use the extrusion method to create a layered mesh. A description of the device with more details is provided in NEGF–Poisson simulations of an FD-SOI field-effect transistor with a quantum dot.

5.5. Setting up the device

5.5.2. Loading the mesh

As usual, the mesh is produced from the appropriate .geo script in the meshes folder, which generates a .msh mesh file and a .xao raw geometry file. Here, we start from an extruded mesh generated by Gmsh. In further steps we will refine the mesh using the adaptive meshing in our Poisson solver. We will then interpolate the resulting potential onto the extruded mesh for the NEGF calculation.

# Set up paths
script_dir = Path(__file__).parent.resolve()
path_mesh = script_dir / "meshes" / "fdsoi_poisson_negf_master.msh"
path_mesh_refined = script_dir / "meshes" / "refined_fdsoi.msh"
path_geo = script_dir / "meshes" / "fdsoi_poisson_negf_master.xao"
path_out = script_dir / "output"

# Load the mesh
scaling = 1e-9
mesh = Mesh(scaling, path_mesh)

Note that we have also defined the path to the folder in which simulation outputs will be saved, path_out.

5.5.3. Defining the device

In this tutorial, we will use two Device objects for the same physical device. The first device (dvc) is used for non-linear Poisson calculations, where the device electrostatics will be solved using adaptive meshing. The second device (dvc_negf) is defined over an extruded mesh along the transport direction to satisfy the NEGF calculation requirements.

We first define variables associated with the gate biases.

# Define the gate parameters
back_gate_bias = -0.2
doping_back_gate = 1e15 * 1e6
binding_energy = 46 * 1e-3 * ct.e
barrier_gate_1 = 0.6
plunger_gate = 0.707
barrier_gate_2 = 0.6
ndoping = 2e20 * 1e6
Ew = mt.Si.Eg/2 + mt.Si.chi # Midgap

Next, we define the Device objects, and then set up the temperature (\(100~\mathrm{mK}\)) and the regions for these devices, in which confined carriers are electrons.

# Define the device object
dvc = Device(mesh, conf_carriers='e')
dvc.set_temperature(0.1)

# Create the regions
dvc.new_region("top_ox", mt.SiO2_ideal)
dvc.new_region("source", mt.Si, ndoping=ndoping)
dvc.new_region("channel", mt.Si)
dvc.new_region("dot", mt.Si)
dvc.new_region("drain", mt.Si, ndoping=ndoping)
dvc.new_region("side_bot_ox", mt.SiO2_ideal)

We define a list of region labels that, taken together, define the channel region which will be used in the NEGF formalism.

Additionally, we define the dot region in which classical charge densities are set to zero. This region will be used later when setting up non-linear Poisson and Schrödinger solvers.

# defining the region forming the quantum dot and the channel.
dvc.set_dot_region("dot")
channel_region_list = ["source", "channel", "dot", "drain"]

# define the boundaries
dvc.new_gate_bnd("barrier_gate_1_bnd", barrier_gate_1, Ew)
dvc.new_gate_bnd("plunger_gate_bnd", plunger_gate, Ew)
dvc.new_gate_bnd("barrier_gate_2_bnd", barrier_gate_2, Ew)
dvc.new_frozen_bnd("back_gate_bnd", back_gate_bias, mt.Si, doping_back_gate,
   "n", binding_energy)

# recreate the device with the extruded mesh for the channel.
dvc_negf = deepcopy(dvc)

In the last line of the above code snippet, we copy the created device with the extruded mesh to be used in the NEGF calculation.

5.6. Simulating the device and results

5.6.1. Configuring the solvers

Next, we configure the parameters for the non-linear Poisson, Schrödinger, NEGF, and lever-arm solvers.

For the non-linear Poisson solver, we use adaptive meshing to facilitate convergence at cryogenic temperature. See the Poisson solver with adaptive meshing and Adaptive-mesh Poisson solver combined with the Schrödinger solver tutorials for a detailed explanation of non-linear Poisson solver parameters.

For the Schrödinger solver, we limit the number of states to four to accelerate the calculation, and to be only slightly larger than the dimension of the basis set that will be employed later in the many-body solver that underpins the master equation calculations of Coulomb peaks and diamonds.

For the NEGF solver, we simply use the default solver parameters. We also define the start and end points of a linecut that crosses the dot inside the channel, which will be used for visualization purposes (e.g. to plot linecuts of the local density of states (LDOS)). For a more detailed NEGF example, refer to the NEGF–Poisson simulations of an FD-SOI field-effect transistor with a quantum dot tutorial.

Finally, the lever arm is an important parameter for the master equation solver. Therefore, we also instantiate a lever-arm matrix Solver object and its parameters.

# setup the Poisson parameters
params_poisson = PoissonSolverParams()
params_poisson.tol = 1e-3 # Convergence threshold (tolerance) for the error
params_poisson.initial_ref_factor = 0.1
params_poisson.final_ref_factor = 0.75
params_poisson.min_nodes = 5e4
params_poisson.max_nodes = 1e5
params_poisson.maxiter_adapt = 100
params_poisson.maxiter = 2000
params_poisson.dot_region = channel_region_list
params_poisson.h_dot = 0.15
params_poisson.mesh_from_file = True
params_poisson.refined_mesh_filename =path_mesh_refined

# Configure the Schrodinger solver
params_schrodinger = SchrodingerSolverParams()
params_schrodinger.tol= 1e-9
params_schrodinger.num_states = 4

# Configure the NEGF-Poisson solver
solver_params = NEGFSolverParams()

# Define parameters pertaining to linecuts over which various quantities will
# be plotted. This line will cross the dot inside the channel.
t_top_ox = 2 * scaling # top oxide thickness
t_si = 6 * scaling # silicon film thickness
begin = (0, np.amin(mesh.glob_nodes[:, 1]), -t_top_ox - t_si/5)
end = (0, np.amax(mesh.glob_nodes[:, 1]), -t_top_ox - t_si/5)
show_figure = True # whether plots are shown

# Solver params for the LeverArmSolver
lam_params = LeverArmMatrixSolverParams()
lam_params.pot_solver_params = params_poisson
lam_params.schrod_solver_params = params_schrodinger

5.6.2. Electrostatics and single-body calculation

In this step, we solve the non-linear Poisson equation in the entire device using the adaptive meshing method. Additionally, we solve the Schrodinger equation in the dot region.

# Instantiate Poisson solver
poisson_slv = PoissonSolver(dvc, solver_params=params_poisson,
   geo_file=path_geo)

# Instantiate Poisson solver
poisson_slv.solve()

# Schrodinger equation calculations.
dvc.set_V_from_phi()

# Creating the dot region subdevice.
submesh = SubMesh(dvc.mesh, ["dot"])
subdevice = SubDevice(dvc, submesh)

# Instantiate Schrodinger solver
schrod_solver = SchrodingerSolver(subdevice, solver_params=params_schrodinger)

# Solve Schrodinger's equation
schrod_solver.solve()

5.6.3. NEGF simulation and boundary conditions

As described in NEGF–Poisson simulations of an FD-SOI field-effect transistor with a quantum dot, all QTCAD NEGF simulations require exactly two labeled surfaces to be assigned a “quantum lead” boundary condition, as may be done using the new_quantum_lead_bnd method of the Device class. These boundaries correspond to the interfaces with the source and drain leads, which, in the NEGF formalism, are assumed to be thermal baths from which charge carriers may be emitted and into which charge carriers may be absorbed.

# define the quantum lead boundary for the NEGF solver
# in the NEGF device.
dvc_negf.new_quantum_lead_bnd("source_bnd", essential=False)
dvc_negf.new_quantum_lead_bnd("drain_bnd", essential=False)

Next, we create the mesh and subdevice for the channel of the FD-SOI device for the NEGF simulation. The NEGF solver requires the entire channel of the device to be provided as the d_negf argument, with the corresponding source and drain regions specified as source_subd and drain_subd.

# Creating a submesh for the channel, source and drain.
submesh = SubMesh(mesh, channel_region_list)
submesh_source = SubMesh(mesh, ["source"])
submesh_drain = SubMesh(mesh, ["drain"])

# Restrict calculation of charge density on a subdevice
d_negf = SubDevice(dvc_negf, submesh)
d_source = SubDevice(dvc_negf, submesh_source)
d_drain = SubDevice(dvc_negf, submesh_drain)

After defining the NEGF device with the layered mesh and specifying the required subdevices, we can interpolate the equilibrium potential values obtained from the non-linear Poisson calculation onto the NEGF device.

# Interpolate the potential for the NEGF device and channel subdevices.
phi_interp = griddata(dvc.mesh.glob_nodes, dvc.phi, mesh.glob_nodes, method='nearest')
phi_interp_n = griddata(dvc.mesh.glob_nodes, dvc.phi, submesh.glob_nodes, method='nearest')
phi_interp_s = griddata(dvc.mesh.glob_nodes, dvc.phi, submesh_source.glob_nodes, method='nearest')
phi_interp_d = griddata(dvc.mesh.glob_nodes, dvc.phi, submesh_drain.glob_nodes, method='nearest')

To simulate near-equilibrium transport, we apply a small, finite source-to-drain bias (Vds). Finally, we initialize the NEGF solver and update the potentials and Hamiltonian of the solver.

# Set up the NEGF solver
Vds = 0.001 # drain-to-source voltage
NEGF_solver = NEGFSolver(d=dvc_negf, d_negf=d_negf, source_subd=d_source,
               drain_subd=d_drain, Vds=Vds, solver_params=solver_params)

# Set the potential from the adaptive mesh and update the Hamiltonian used in NEGF
NEGF_solver.d.set_potential(phi_interp)
NEGF_solver.d_negf.set_potential(phi_interp_n)
NEGF_solver.source_subd.set_potential(phi_interp_s)
NEGF_solver.drain_subd.set_potential(phi_interp_d)

# Update the solver and the Hamiltonian
NEGF_solver.update_hamiltonian()

After updating the potentials and the Hamiltonian, we can calculate the current.

# Compute the current
path_current = str(path_out / ("fdsoi_poisson_negf_master_current.png"))
current = NEGF_solver.get_current_adaptive(log_scale=False,
   show_figure=show_figure, path=path_current,
   title="Spectral current near equilibrium")
print("Current: "+str(current)+" [A]")

Executing the above lines of code results in the following output.

================================================================================
NEGF-POISSON SOLVER INITIALIZED
Current: 6.13449069876267e-09 [A]

Additionally, in the above lines of code, we plot the spectral current as a function of energy.

Fig. 5.6.1 Transmission function in the sequential tunneling regime.

Inspecting the figure, we notice a sharp peak which corresponds to a quasi-bound state. This is expected, as we want to operate in the sequential tunneling regime.

To better illustrate the physics of the tunneling and transport phenomena we can plot the LDOS (see Eq. (11.12)). We also plot the LDOS with narrower bounds on the energy axis to better visualize the LDOS between the source and drain Fermi levels.

# compute LDOS
path_ldos = str(path_out / ("fdsoi_poisson_negf_master_ldos.png"))
NEGF_solver.plot_ldos(begin=begin, end=end, title="Local density of states",
   show_figure=show_figure, path=path_ldos, energies = np.linspace(-0.1*ct.e, 0.03*ct.e, 50))

# compute LDOS (zoomed-in)
path_ldos = str(path_out / ("fdsoi_poisson_negf_master_ldos_zoom.png"))
NEGF_solver.plot_ldos(begin=begin, end=end, title="Local density of states",
   show_figure=show_figure, path=path_ldos, energies = np.linspace(-0.02*ct.e, 0.01*ct.e, 50))

This results in the following figures.

Fig. 5.6.2 LDOS plot in the sequential tunneling regime.

Fig. 5.6.3 LDOS plot in the sequential tunneling regime zoomed-in on the quasi-bound state.

We notice that the first quasi-bound state localized inside the dot (which corresponds to the quantum dot ground state) lies between the source and drain Fermi levels. Therefore, the current value through the channel corresponds to the amplitude of the Coulomb peak associated with electrons tunneling through the ground state of the quantum dot. Essentially, the bias configuration chosen here leads to tunneling from the source (left lead) to the dot and from the dot to the right lead (drain).

5.6.4. Solving the master equation

After finding the sequential-tunneling-regime current from the NEGF simulation, we can calculate the Coulomb-blockade-regime current by solving the master equation within the featureless approximation.

To accelerate the calculation, we calculate the lever arm of the plunger gate on the quantum dot.

# Instantiate lever arm matrix solver
slv = LeverArmMatrixSolver(dvc, ["plunger_gate_bnd"], [plunger_gate], dot_region="dot",
   solver_params=lam_params)

# Calculate the lever arm matrix
bias_increment = 1e-3
lever_arm_matrix = slv.solve(bias_increment=bias_increment)

After finding the lever arm of the plunger gate, we configure the many-body solver and instantiate a Junction object to obtain the current by solving the master equation. In this tutorial, we only consider two spin-degenerate single-body basis states as we set num_states = 2.

# Configure the many-body solver for the junction class.
params_mb = ManyBodySolverParams()
params_mb.n_degen = 2
params_mb.num_states = 2   # A larger value can be used for more accurate results
params_mb.alpha = lever_arm_matrix[0,0]

# Initializing the Junction object used to calculate transport
jc = Junction(subdevice, many_body_solver_params=params_mb)

In this step, we apply the same source–drain bias used in the NEGF calculation and we retrieve the Coulomb peak positions from the coulomb_peak_pos attribute of the Junction object.

# Set a near-zero source-drain bias (same as the Vds used in NEGF solver)
jc.setVs(Vds)
jc.setVd(0)
# printing the Coulomb peak positions
print('Coulomb peak positions:')
print(jc.coulomb_peak_pos)

Since the first localized bound state in the LDOS plot is aligned with the source and drain Fermi levels, we expect the first Coulomb peak position to be near zero as our chosen reference bias configuration corresponds to the first Coulomb peak. This leads to the following output.

Coulomb peak positions:
[-0.02268238  0.05472169  0.14125148  0.20845242]

We notice that the first peak position is quite close to zero. However, this small offset can be explained by the following points:

  • Since the two calculations are based on different methods, numerical errors are expected to cause deviations.

  • Tunnel coupling of the channel with the leads may also cause small shifts in the peak positions that are not captured in the many-body solver that underpins the master equation technique. Since this solver is defined in the quantum dot region exclusively, the eigenstates do not account for wavefunction penetration into the leads, which is considered in the NEGF formalism.

In the next step, we solve the master equation to find the current through the channel over a range of bias values for the plunger gate potential (using the plunger-gate lever arm). The range of the bias sweep is determined from the Coulomb peak positions.

gate_bias_sweep = np.linspace(-0.1, 0.3, 400)
vec_Il = []
bartitle = "Calculating Coulomb peaks."
progress_bar = Bar(bartitle, max = len(gate_bias_sweep))
for idx, val in enumerate(gate_bias_sweep):
   jc.setVg(val)
   Il,Ir,prob = seq_tunnel_curr(jc)    # transport calculation
   vec_Il += [Il]
   progress_bar.next()
progress_bar.finish()

vec_Il = np.array(vec_Il)
np.savetxt(path_out/'PNEGF_current_scale.txt',vec_Il)

5.6.5. Coulomb peak results

Using the following lines of code, we can plot the current values obtained using the master equation. Since the broadening function used to obtain these values within the featureless approximation is an arbitrary number (set by default to \(1~\mathrm{Hz}\)) the units of the current are also arbitrary.

peaks, _  = find_peaks(vec_Il, 2e-20)
fig, ax = plt.subplots()
ax.plot(gate_bias_sweep, vec_Il, '-', label='Master equation')
ax.plot(gate_bias_sweep[peaks], vec_Il[peaks],'x', label= 'Coulomb peak maxima')
ax.set_ylabel('Current (arbitrary units)')
ax.set_xlabel('Gate bias (V)')
ax.legend()
fig.savefig(path_out/'PNEGF_current_featureless.png')
plt.show()

Fig. 5.6.4 Coulomb peaks in arbitrary units

Since we have previously obtained the current value corresponding to the maximum of the first peak using the NEGF formalism, we can scale the currents obtained in the master equation method to match the maximum of the first peak with the NEGF current.

# scaling the values
scaling_value =  current/vec_Il[peaks[0]]
vec_Il = scaling_value*vec_Il

We may then plot the scaled currents obtained in the master equation method with the following code.

fig, ax = plt.subplots()
ax.plot(gate_bias_sweep, vec_Il, '-')
ax.plot(gate_bias_sweep[peaks], vec_Il[peaks],'x')
ax.set_ylabel('Current (A)')
ax.set_xlabel('Gate bias (V)')
fig.savefig(path_out/'PNEGF_Current_scale.png')
plt.show()

This results in the following plot.

Fig. 5.6.5 Coulomb peaks in SI units.

Similar to the previous section, we can calculate the charge stability diagram using the master equation and rescale the currents to obtain physically-meaningful values.

# Calculating the charge stability diagram
source_drain_sweep = np.linspace(-0.1, 0.1, 150)
gate_bias_sweep = np.linspace(-0.1, 0.3, 150)
grid_Il = np.zeros((150, 150))
bartitle = "Calculating charge-stability diagram."
progress_bar = Bar(bartitle, max = len(gate_bias_sweep)*len(source_drain_sweep))
for idx, val_p in enumerate(gate_bias_sweep):
   jc.setVg(val_p)
   for idy, val in enumerate(source_drain_sweep):
      jc.setVs(val)
      Il,Ir,prob = seq_tunnel_curr(jc)    # transport calculation
      grid_Il[idx, idy] = Il
      progress_bar.next()
progress_bar.finish()

# Saving the data
np.savetxt(path_out/'PNEGF_CSD_scale.txt',grid_Il)

# Plotting the charge stability diagram
fig, ax = plt.subplots()
grid_Il = scaling_value*grid_Il
im = ax.imshow(np.transpose(grid_Il), interpolation='bilinear',
   cmap='bwr', extent=[-0.1, 0.3, -0.1, 0.1], aspect=2)
ax.set_ylabel(r'Drain--source bias $V_{DS}$ (V)')
ax.set_xlabel('Gate bias (V)')
cbar = fig.colorbar(im)
cbar.ax.set_ylabel('Current (A)')
fig.savefig(path_out/'PNEGF_Charge_diagram.png')
plt.show()

This results in the following charge stability diagram.

Fig. 5.6.6 Charge stability diagram in SI units.

5.7. Full code

__copyright__ = "Copyright 2022-2024, Nanoacademic Technologies Inc."

import numpy as np
from pathlib import Path
from progress.bar import ChargingBar as Bar
from matplotlib import pyplot as plt
from copy import deepcopy
from scipy.interpolate import griddata
from scipy.signal import find_peaks
from qtcad.device import constants as ct
from qtcad.device import materials as mt
from qtcad.device.mesh3d import Mesh, SubMesh
from qtcad.device import Device, SubDevice
from qtcad.device.poisson import Solver as PoissonSolver
from qtcad.device.poisson import SolverParams as PoissonSolverParams
from qtcad.device.schrodinger import Solver as SchrodingerSolver
from qtcad.device.schrodinger import SolverParams as SchrodingerSolverParams
from qtcad.transport.negf_poisson import Solver as NEGFSolver
from qtcad.transport.negf_poisson import SolverParams as NEGFSolverParams
from qtcad.device.leverarm_matrix import Solver as LeverArmMatrixSolver
from qtcad.device.leverarm_matrix import SolverParams as LeverArmMatrixSolverParams
from qtcad.device.many_body import SolverParams as ManyBodySolverParams
from qtcad.transport.junction import Junction
from qtcad.transport.mastereq import seq_tunnel_curr

# Set up paths
script_dir = Path(__file__).parent.resolve()
path_mesh = script_dir / "meshes" / "fdsoi_poisson_negf_master.msh"
path_mesh_refined = script_dir / "meshes" / "refined_fdsoi.msh"
path_geo = script_dir / "meshes" / "fdsoi_poisson_negf_master.xao"
path_out = script_dir / "output"

# Load the mesh
scaling = 1e-9
mesh = Mesh(scaling, path_mesh)

# Define the gate parameters
back_gate_bias = -0.2
doping_back_gate = 1e15 * 1e6
binding_energy = 46 * 1e-3 * ct.e
barrier_gate_1 = 0.6
plunger_gate = 0.707
barrier_gate_2 = 0.6
ndoping = 2e20 * 1e6
Ew = mt.Si.Eg/2 + mt.Si.chi # Midgap

# Define the device object
dvc = Device(mesh, conf_carriers='e')
dvc.set_temperature(0.1)

# Create the regions
dvc.new_region("top_ox", mt.SiO2_ideal)
dvc.new_region("source", mt.Si, ndoping=ndoping)
dvc.new_region("channel", mt.Si)
dvc.new_region("dot", mt.Si)
dvc.new_region("drain", mt.Si, ndoping=ndoping)
dvc.new_region("side_bot_ox", mt.SiO2_ideal)

# defining the region forming the quantum dot and the channel.
dvc.set_dot_region("dot")
channel_region_list = ["source", "channel", "dot", "drain"]

# define the boundaries
dvc.new_gate_bnd("barrier_gate_1_bnd", barrier_gate_1, Ew)
dvc.new_gate_bnd("plunger_gate_bnd", plunger_gate, Ew)
dvc.new_gate_bnd("barrier_gate_2_bnd", barrier_gate_2, Ew)
dvc.new_frozen_bnd("back_gate_bnd", back_gate_bias, mt.Si, doping_back_gate, 
   "n", binding_energy)

# recreate the device with the extruded mesh for the channel.
dvc_negf = deepcopy(dvc)

# setup the Poisson parameters
params_poisson = PoissonSolverParams()
params_poisson.tol = 1e-3 # Convergence threshold (tolerance) for the error
params_poisson.initial_ref_factor = 0.1
params_poisson.final_ref_factor = 0.75
params_poisson.min_nodes = 5e4
params_poisson.max_nodes = 1e5
params_poisson.maxiter_adapt = 100
params_poisson.maxiter = 2000
params_poisson.dot_region = channel_region_list
params_poisson.h_dot = 0.15
params_poisson.mesh_from_file = True
params_poisson.refined_mesh_filename =path_mesh_refined

# Configure the Schrodinger solver
params_schrodinger = SchrodingerSolverParams()
params_schrodinger.tol= 1e-9
params_schrodinger.num_states = 4

# Configure the NEGF-Poisson solver
solver_params = NEGFSolverParams()

# Define parameters pertaining to linecuts over which various quantities will
# be plotted. This line will cross the dot inside the channel.
t_top_ox = 2 * scaling # top oxide thickness
t_si = 6 * scaling # silicon film thickness
begin = (0, np.amin(mesh.glob_nodes[:, 1]), -t_top_ox - t_si/5)
end = (0, np.amax(mesh.glob_nodes[:, 1]), -t_top_ox - t_si/5)
show_figure = True # whether plots are shown

# Solver params for the LeverArmSolver
lam_params = LeverArmMatrixSolverParams()
lam_params.pot_solver_params = params_poisson
lam_params.schrod_solver_params = params_schrodinger

# Instantiate Poisson solver
poisson_slv = PoissonSolver(dvc, solver_params=params_poisson, 
   geo_file=path_geo)

# Instantiate Poisson solver
poisson_slv.solve()

# Schrodinger equation calculations.
dvc.set_V_from_phi()

# Creating the dot region subdevice.
submesh = SubMesh(dvc.mesh, ["dot"])
subdevice = SubDevice(dvc, submesh)

# Instantiate Schrodinger solver
schrod_solver = SchrodingerSolver(subdevice, solver_params=params_schrodinger)

# Solve Schrodinger's equation
schrod_solver.solve()

# define the quantum lead boundary for the NEGF solver
# in the NEGF device.
dvc_negf.new_quantum_lead_bnd("source_bnd", essential=False)
dvc_negf.new_quantum_lead_bnd("drain_bnd", essential=False)

# Creating a submesh for the channel, source and drain.
submesh = SubMesh(mesh, channel_region_list)
submesh_source = SubMesh(mesh, ["source"])
submesh_drain = SubMesh(mesh, ["drain"])

# Restrict calculation of charge density on a subdevice
d_negf = SubDevice(dvc_negf, submesh)
d_source = SubDevice(dvc_negf, submesh_source)
d_drain = SubDevice(dvc_negf, submesh_drain)

# Interpolate the potential for the NEGF device and channel subdevices.
phi_interp = griddata(dvc.mesh.glob_nodes, dvc.phi, mesh.glob_nodes, method='nearest')
phi_interp_n = griddata(dvc.mesh.glob_nodes, dvc.phi, submesh.glob_nodes, method='nearest')
phi_interp_s = griddata(dvc.mesh.glob_nodes, dvc.phi, submesh_source.glob_nodes, method='nearest')
phi_interp_d = griddata(dvc.mesh.glob_nodes, dvc.phi, submesh_drain.glob_nodes, method='nearest')

# Set up the NEGF solver
Vds = 0.001 # drain-to-source voltage
NEGF_solver = NEGFSolver(d=dvc_negf, d_negf=d_negf, source_subd=d_source, 
               drain_subd=d_drain, Vds=Vds, solver_params=solver_params)

# Set the potential from the adaptive mesh and update the Hamiltonian used in NEGF
NEGF_solver.d.set_potential(phi_interp)
NEGF_solver.d_negf.set_potential(phi_interp_n)
NEGF_solver.source_subd.set_potential(phi_interp_s)
NEGF_solver.drain_subd.set_potential(phi_interp_d)

# Update the solver and the Hamiltonian
NEGF_solver.update_hamiltonian()

# Compute the current
path_current = str(path_out / ("fdsoi_poisson_negf_master_current.png"))
current = NEGF_solver.get_current_adaptive(log_scale=False,
   show_figure=show_figure, path=path_current,
   title="Spectral current near equilibrium")
print("Current: "+str(current)+" [A]")

# compute LDOS
path_ldos = str(path_out / ("fdsoi_poisson_negf_master_ldos.png"))
NEGF_solver.plot_ldos(begin=begin, end=end, title="Local density of states",
   show_figure=show_figure, path=path_ldos, energies = np.linspace(-0.1*ct.e, 0.03*ct.e, 50))

# compute LDOS (zoomed-in)
path_ldos = str(path_out / ("fdsoi_poisson_negf_master_ldos_zoom.png"))
NEGF_solver.plot_ldos(begin=begin, end=end, title="Local density of states",
   show_figure=show_figure, path=path_ldos, energies = np.linspace(-0.02*ct.e, 0.01*ct.e, 50))

# Instantiate lever arm matrix solver
slv = LeverArmMatrixSolver(dvc, ["plunger_gate_bnd"], [plunger_gate], dot_region="dot",
   solver_params=lam_params)

# Calculate the lever arm matrix
bias_increment = 1e-3
lever_arm_matrix = slv.solve(bias_increment=bias_increment)

# Configure the many-body solver for the junction class.
params_mb = ManyBodySolverParams()
params_mb.n_degen = 2
params_mb.num_states = 2   # A larger value can be used for more accurate results
params_mb.alpha = lever_arm_matrix[0,0]

# Initializing the Junction object used to calculate transport
jc = Junction(subdevice, many_body_solver_params=params_mb)

# Set a near-zero source-drain bias (same as the Vds used in NEGF solver)
jc.setVs(Vds)
jc.setVd(0)
# printing the Coulomb peak positions
print('Coulomb peak positions:')
print(jc.coulomb_peak_pos)

gate_bias_sweep = np.linspace(-0.1, 0.3, 400)
vec_Il = []
bartitle = "Calculating Coulomb peaks."
progress_bar = Bar(bartitle, max = len(gate_bias_sweep))
for idx, val in enumerate(gate_bias_sweep):
   jc.setVg(val)
   Il,Ir,prob = seq_tunnel_curr(jc)    # transport calculation
   vec_Il += [Il]
   progress_bar.next()
progress_bar.finish()

vec_Il = np.array(vec_Il)
np.savetxt(path_out/'PNEGF_current_scale.txt',vec_Il)

peaks, _  = find_peaks(vec_Il, 2e-20)
fig, ax = plt.subplots()
ax.plot(gate_bias_sweep, vec_Il, '-', label='Master equation')
ax.plot(gate_bias_sweep[peaks], vec_Il[peaks],'x', label= 'Coulomb peak maxima')
ax.set_ylabel('Current (arbitrary units)')
ax.set_xlabel('Gate bias (V)')
ax.legend()
fig.savefig(path_out/'PNEGF_current_featureless.png')
plt.show()

# scaling the values
scaling_value =  current/vec_Il[peaks[0]]
vec_Il = scaling_value*vec_Il

fig, ax = plt.subplots()
ax.plot(gate_bias_sweep, vec_Il, '-')
ax.plot(gate_bias_sweep[peaks], vec_Il[peaks],'x')
ax.set_ylabel('Current (A)')
ax.set_xlabel('Gate bias (V)')
fig.savefig(path_out/'PNEGF_Current_scale.png')
plt.show()

# Calculating the charge stability diagram
source_drain_sweep = np.linspace(-0.1, 0.1, 150)
gate_bias_sweep = np.linspace(-0.1, 0.3, 150)
grid_Il = np.zeros((150, 150))
bartitle = "Calculating charge-stability diagram."
progress_bar = Bar(bartitle, max = len(gate_bias_sweep)*len(source_drain_sweep))
for idx, val_p in enumerate(gate_bias_sweep):
   jc.setVg(val_p)
   for idy, val in enumerate(source_drain_sweep):
      jc.setVs(val)
      Il,Ir,prob = seq_tunnel_curr(jc)    # transport calculation
      grid_Il[idx, idy] = Il
      progress_bar.next()
progress_bar.finish()

# Saving the data
np.savetxt(path_out/'PNEGF_CSD_scale.txt',grid_Il)

# Plotting the charge stability diagram
fig, ax = plt.subplots()
grid_Il = scaling_value*grid_Il
im = ax.imshow(np.transpose(grid_Il), interpolation='bilinear', 
   cmap='bwr', extent=[-0.1, 0.3, -0.1, 0.1], aspect=2)
ax.set_ylabel(r'Drain--source bias $V_{DS}$ (V)')
ax.set_xlabel('Gate bias (V)')
cbar = fig.colorbar(im)
cbar.ax.set_ylabel('Current (A)')
fig.savefig(path_out/'PNEGF_Charge_diagram.png')
plt.show()