4. Computing the chemical potential of an SET device

4.1. Requirements

4.1.1. Software components

  • QTCAD

4.1.2. Python script

  • qtcad/examples/practical_application/FDSOI/chemical_potential.py

4.1.3. References

4.2. Briefing

In this section of the tutorial, we will demonstrate how to compute the chemical potential, \(\mu\), of the single-electron transistor (SET) as a function of the gate voltage applied to the associated gate, in this case, "plunger_gate_1_bnd".

In principle, to properly compute the chemical potential, one would need to solve the full many-body Schrödinger equation. However, this is computationally very expensive and not feasible for large systems. Instead, we will use a simplified approach to account for many-body effects by solving the Schrödinger–Poisson equation, which self-consistently couples the single-particle Schrödinger equation with the Poisson equation. This approach captures the essential physics of electron-electron interactions in a mean-field manner, making it computationally more tractable.

Note

While the Schrödinger–-Poisson approach (outlined in Schrödinger–Poisson solver) incorporates many-body effects in a mean-field manner, it has important limitations. In particular, it does not account for exchange and correlation contributions to the energy. Consequently, results obtained with this method should be interpreted with caution in regimes where strong correlations are expected to be significant.

Furthermore, unlike approaches such as Hartree–-Fock [BSA00], the Schrödinger–-Poisson method described in Schrödinger–Poisson solver does not eliminate electron self-interaction, which can introduce inaccuracies in the computed chemical potentials. We expect the relative impact of self-interaction errors to diminish as the number of electrons increases, making this approach more reliable for systems with a large number of electrons.

4.4. Solver paramaters

The main goal of this tutorial series is to simulate the behavior of an SET using the Schrödinger–Poisson solver. Because most of the steps of the tutorial involve solving the Schrödinger–Poisson equations, we define here the default SolverParams parameters that will be used in the calculations.

# Solver parameters
# (Poisson)
p_params                            = PoissonSolverParams()
p_params.tol                        = 1e-8
# (Schrödinger)
schrod_params                       = SchrodSolverParams()
schrod_params.tol                   = 1e-8
schrod_params.num_states            = 10
# (Schrödinger-Poisson)
sp_params                           = SolverParams()
sp_params.tol                       = 1e-5
sp_params.maxiter                   = 2000
sp_params.bound_state_charges_only  = True
sp_params.sc_method                 = 'underrelax'
sp_params.initialization            = False
sp_params.mixing_algo               = "adaptive_linear"
sp_params.mixing_param              = 0.1
sp_params.poisson_solver_params     = p_params
sp_params.schrod_solver_params      = schrod_params

We start by setting up the parameters for the Poisson solver (PoissonSolverParams) and the Schrödinger solver (SchrodSolverParams). Then we create an instance of the Schrödinger–Poisson SolverParams class that contains the parameters for the Schrödinger–Poisson Solver. These parameters include:

  • tol — Convergence tolerance for the self-consistent loop. Here the tolerance is set to 1e-5 eV, meaning that the solver will consider the solution converged when the maximal relative change in the potential between successive iterations is less than this value. The reason we require such a strict tolerance is that we will later compute many-body energies, chemical potentials, and addition energies which are expressed as sums or difference of single-particle eigenenergies. The typical values we will be interested in are on the order of one to ten meV. Because we will be interested in energies of \(N \sim 12\) electrons, the error on the quantities of interest can be on the order of 10 times larger than the we error on the single-particle eigenenergies. Therefore, we need an accuracy of each single-particle eigenenergy to be significantly better than 1 meV to ensure that the computed energies are accurate.

  • maxiter — Maximum number of iterations for the self-consistent loop.

  • bound_state_charges_only — Whether or not to consider only bound-state charges. Here it is set to True, specifying that the only charges contributing to the Poisson equation are those from bound states (electrons occupying eigenstates of the Schrödinger equation). If False, classical charges from conduction and valence bands would also be included.

  • sc_method — Self-consistent method to be used. Here it is set to 'underrelax'. The other option is 'p-c' for predictor–corrector. While predictor–corrector has been shown to be more efficient when simulating systems containing electron gases [TGPR97] originating from the filling of parabolic bands, here, all such charge is neglected (see description of bound_state_charges_only parameter, above).

  • initialization — Whether to generate an initial guess for the electric potential phi. Here, it is set to False since we will provide our own initial guess, either by first solving the linear Poisson equation (see Simulating electrostatics using the linear Poisson solver) or by using a previously computed solution.

  • mixing_algo — Mixing algorithm to be used. Here, we use the "adaptive_linear" algorithm which adaptively adjusts the mixing parameter for a linear-based mixing scheme, based on the convergence behavior.

  • mixing_param — Mixing parameter. In the case of the "adaptive_linear" algorithm, this parameter serves as the initial mixing parameter.

  • poisson_solver_params — The Poisson solver parameters defined earlier.

  • schrod_solver_params — The Schrödinger solver parameters defined earlier.

4.5. Function to compute the chemical potential

Next, we define the main content of the corresponding Python file: a function, compute_chemical_potential, which computes the chemical potential as a function of the gate voltage applied to the SET gate, "plunger_gate_1_bnd".

4.5.1. Function description

The first step is to define the function signature and docstrings.

def compute_chemical_potential(V_SET: np.ndarray, Nmin: int, Nmax: int,
    d: Device, set_device: SubDevice, out_path: pathlib.Path=script_dir) -> np.ndarray:
    """Compute the chemical potential for different number of particles and
    different voltage configurations. This function also saves the potential
    and wavefunctions to hdf5 files, and the energies and population factors
    to text files. All of the output files are saved in the specified output
    directory.

    Args:
        V_SET (np.ndarray): Array of SET gate voltages.
        Nmin (int): Minimum number of electrons.
        Nmax (int): Maximum number of electrons.
        d (Device): Device object representing the entire FD-SOI device.
        set_device (SubDevice): SubDevice representing the SET region of the
            FD-SOI device.
        out_path (pathlib.Path): Output directory path.

    Returns:
        np.ndarray: 2D array of chemical potentials for each voltage
        configuration. The first index runs over the applied voltages to the
        SET gate, and the second index runs over the number of particles.

    Note:
        The chemical potential is calculated as the difference in many-body
        ground-state energies when adding an electron to the system. The
        chemical potentials are computed for the electron numbers in the range
        [Nmin+1, Nmax].
    """

As is explained in the docstring, the function takes as input an array of gate voltages, the minimum and maximum number of electrons to consider, the Device object representing the entire FD-SOI device, the SubDevice object representing the SET region, and the output directory path. It returns a 2D array of chemical potentials for each voltage configuration and number of particles.

As mentioned above, in this tutorial, we will compute the chemical potential when the number of electrons in the SET is around \(N=12\), for different values of the "plunger_gate_1_bnd" voltage. The chemical potential \(\mu(N)\) is defined as the energy required to add the \(N\)-th electron to the system, i.e., \(\mu(N) = U(N) - U(N-1)\), where \(U(N)\) is the ground-state energy of the system with \(N\) electrons [HKP+07]. Therefore, for each value of the gate voltage, we will compute the total energy of the system by converging a Schrödinger–Poisson calculation for \(N-1 = 12\), \(N = 13\), and \(N+1=14\) electrons to obtain the chemical potentials \(\mu(N=13)\) and \(\mu(N=14)\). Users are free to modify the range of electron numbers and gate voltages as needed.

4.5.2. Preliminaries

Inside the function, the first step is to define the paths to the output files that will be used to save the results of the simulations

# output files
# (results)
phi_out_file         = str(out_path / "phi_vset{vset:.2f}_N{n}.hdf5")
wf_hdf5_file         = str(out_path / "wf_vset{vset:.2f}_N{n}.hdf5")
mu_file              = out_path / "mu.txt"
pop_file             = out_path / "population.txt"
E_file               = out_path / "energies.txt"

# (visualization)
wf_vtu_file          = str(out_path / "wf_vset{vset:.2f}_N{n}.vtu")
rho_slice_file       = str(out_path / "rho_vset{vset:.2f}_N{n}_slice.png")
wf_slice_file        = str(out_path / "wf_vset{vset:.2f}_N{n}_slice.png")

Then, the outputs of intermediary steps and of the entire calculation are initialized.

# Range of number of particles
N = range(Nmin, Nmax + 1)
# Many-body energies
mb_E = np.zeros((len(V_SET), len(N)))
# Chemical potential
Mu = np.zeros((len(V_SET), len(N) - 1))

We start by defining the range of number of particles to consider, from Nmin to Nmax. We then initialize a 2D array, mb_E, to store the many-body ground-state energies for each value of the gate voltage and number of particles. Finally, we initialize a 2D array, Mu, to store the chemical potentials for each value of the gate voltage and number of particles.

4.5.3. Solving the Schrödinger–Poisson equation

Next, we loop over the values of the gate voltage, print some information about which calculation is currently being performed, and set the potential on the SET gate, "plunger_gate_1_bnd", accordingly

for j, vset in enumerate(V_SET):
    # Print info
    print("\n")
    print("==================================================")
    print(f"Setting V_SET = {vset} V")
    print("==================================================")
    print("\n")

    # Set the potential in the SET region
    d.set_applied_potential("plunger_gate_1_bnd", vset)

For each of these applied voltages, we then loop over the number of particles and perform a Schrödinger–Poisson calculation to obtain the ground-state energy of the system

for i, n in enumerate(N):
    # Print info
    print("--------------------------------------------------")
    print(f"Solving for N = {n} electrons...")
    print("--------------------------------------------------")
    print("\n")

    # Solve Schrödinger-Poisson in the SET region
    sp_solver = Solver(d, subdevice=set_device, solver_params=sp_params)
    sp_solver.solve(N=n)
    print("Energies and population factors of SET SubDevice:")
    set_device.print_energies()

The argument N passed to the solve method specifies the number of electrons to be considered in the Schrödinger–Poisson calculation. To enforce the desired number of electrons, the solver places electrons in the single-particle energy levels according to the Fermi-Dirac distribution with the Fermi level adjusted such that the system contains the desired number of electrons.

4.5.4. Saving results

After solving the Schrödinger–Poisson equation, we save the population factors and the single-particle energies to output files

# Save energies and population factors
with open(pop_file, "a") as f:
    out = [vset, n] + list(set_device.population_factors)
    np.savetxt(f, np.array(out)[np.newaxis, :], fmt="%.6f")
with open(E_file, "a") as f:
    out = [vset, n] + list(set_device.energies / ct.e)
    np.savetxt(f, np.array(out)[np.newaxis, :], fmt="%.8f")

The population factors indicate the occupation of each single-particle energy level in the SET region. Each of the factors is a real number between 0 and 1, where a value of 1 indicates that the corresponding energy level is fully occupied, while a value of 0 indicates that the level is empty. It is important to note that the population factors do not include the degeneracy. In other words, each population factor should be multiplied by the degeneracy of the corresponding energy level to obtain the actual number of electrons occupying that level. For example, as mentioned above, in silicon, each energy level has a valley degeneracy of 2 and a spin degeneracy of 2, resulting in a total degeneracy of g=4 (defined above).

Next, we save the potential to an .hdf5 file for later use, the wavefunctions to a .vtu file for visualization, and save slices of the ground-state density and total electron density for quick visualization purposes

# Save phi, wavefunctions, and density
io.save(phi_out_file.format(vset=vset, n=n), {"phi": d.phi})
wf_dict = {f"WF_{i}": np.abs(set_device.eigenfunctions[:, i])**2 for i in range(set_device.eigenfunctions.shape[1])}
io.save(wf_vtu_file.format(vset=vset, n=n), wf_dict, set_device.mesh)

plot_slice(set_device.mesh, np.abs(set_device.eigenfunctions[:, 0])**2,
        normal=normal, origin=origin,
        cb_axis_label="Ground-state density [$e$ / m$^3$]",
        path=wf_slice_file.format(vset=vset, n=n), show_figure=False)
plot_slice(d.mesh, d.rho / ct.e, normal=normal, origin=origin,
        cb_axis_label="$\\rho$ [$e$ / m$^3$]",
        path=rho_slice_file.format(vset=vset, n=n), show_figure=False)

4.5.5. Computing chemical potentials

The final step is to compute the many-body ground-state energy and the chemical potential. The many-body ground-state energy is computed as the sum of the single-particle energies weighted by their respective population factors.

# Compute many-body ground-state energy
mb_E[j, i] = np.sum(g * set_device.population_factors * set_device.energies)
print("  Many-body ground-state energy (eV): "f"{mb_E[j, i] / ct.e:.10f}")

The chemical potential is then computed as the difference in many-body ground-state energies when adding an electron to the system [HKP+07].

        if i > 0:
            # Compute chemical potential
            Mu[j, i-1] = (mb_E[j, i] - mb_E[j, i-1])
            print(f"  Chemical potentials (eV): {Mu[j, i-1] / ct.e}")
            with open(mu_file, "a") as f:
                np.savetxt(f, np.array([[vset, n, Mu[j, i-1] / ct.e]]), fmt="%.10f")

return Mu

4.6. Full code

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

import pathlib
import numpy as np
from qtcad.device import constants as ct
from qtcad.device import materials as mt
from qtcad.device import Device, SubDevice
from qtcad.device.schrodinger_poisson import Solver 
from qtcad.device.schrodinger_poisson import SolverParams
from qtcad.device.schrodinger import SolverParams as SchrodSolverParams
from qtcad.device.poisson import SolverParams as PoissonSolverParams
from qtcad.device import io
from qtcad.device.analysis import plot_slice
from double_dot_fdsoi import normal, origin

# Script directory
script_dir = pathlib.Path(__file__).parent.resolve()

# Degeneracy factor for silicon
g = int(mt.Si.gqc)    

# Solver parameters
# (Poisson)
p_params                            = PoissonSolverParams()
p_params.tol                        = 1e-8
# (Schrödinger)
schrod_params                       = SchrodSolverParams()
schrod_params.tol                   = 1e-8
schrod_params.num_states            = 10
# (Schrödinger-Poisson)
sp_params                           = SolverParams()
sp_params.tol                       = 1e-5
sp_params.maxiter                   = 2000
sp_params.bound_state_charges_only  = True
sp_params.sc_method                 = 'underrelax'
sp_params.initialization            = False
sp_params.mixing_algo               = "adaptive_linear"
sp_params.mixing_param              = 0.1
sp_params.poisson_solver_params     = p_params
sp_params.schrod_solver_params      = schrod_params

def compute_chemical_potential(V_SET: np.ndarray, Nmin: int, Nmax: int, 
    d: Device, set_device: SubDevice, out_path: pathlib.Path=script_dir) -> np.ndarray:
    """Compute the chemical potential for different number of particles and
    different voltage configurations. This function also saves the potential
    and wavefunctions to hdf5 files, and the energies and population factors
    to text files. All of the output files are saved in the specified output
    directory.

    Args:
        V_SET (np.ndarray): Array of SET gate voltages.
        Nmin (int): Minimum number of electrons.
        Nmax (int): Maximum number of electrons.
        d (Device): Device object representing the entire FD-SOI device.
        set_device (SubDevice): SubDevice representing the SET region of the 
            FD-SOI device.
        out_path (pathlib.Path): Output directory path.

    Returns:
        np.ndarray: 2D array of chemical potentials for each voltage 
        configuration. The first index runs over the applied voltages to the
        SET gate, and the second index runs over the number of particles.

    Note:
        The chemical potential is calculated as the difference in many-body 
        ground-state energies when adding an electron to the system. The 
        chemical potentials are computed for the electron numbers in the range
        [Nmin+1, Nmax].
    """

    # output files
    # (results)
    phi_out_file         = str(out_path / "phi_vset{vset:.2f}_N{n}.hdf5")
    wf_hdf5_file         = str(out_path / "wf_vset{vset:.2f}_N{n}.hdf5")
    mu_file              = out_path / "mu.txt"
    pop_file             = out_path / "population.txt"
    E_file               = out_path / "energies.txt"

    # (visualization)
    wf_vtu_file          = str(out_path / "wf_vset{vset:.2f}_N{n}.vtu")
    rho_slice_file       = str(out_path / "rho_vset{vset:.2f}_N{n}_slice.png")
    wf_slice_file        = str(out_path / "wf_vset{vset:.2f}_N{n}_slice.png")

    # Range of number of particles
    N = range(Nmin, Nmax + 1)
    # Many-body energies
    mb_E = np.zeros((len(V_SET), len(N)))
    # Chemical potential
    Mu = np.zeros((len(V_SET), len(N) - 1))

    for j, vset in enumerate(V_SET):
        # Print info
        print("\n")
        print("==================================================")
        print(f"Setting V_SET = {vset} V")
        print("==================================================")
        print("\n")

        # Set the potential in the SET region
        d.set_applied_potential("plunger_gate_1_bnd", vset)

        for i, n in enumerate(N):
            # Print info
            print("--------------------------------------------------")
            print(f"Solving for N = {n} electrons...")
            print("--------------------------------------------------")
            print("\n")

            # Solve Schrödinger-Poisson in the SET region
            sp_solver = Solver(d, subdevice=set_device, solver_params=sp_params)
            sp_solver.solve(N=n)
            print("Energies and population factors of SET SubDevice:")
            set_device.print_energies()

            # Save energies and population factors
            with open(pop_file, "a") as f:
                out = [vset, n] + list(set_device.population_factors)
                np.savetxt(f, np.array(out)[np.newaxis, :], fmt="%.6f")
            with open(E_file, "a") as f:
                out = [vset, n] + list(set_device.energies / ct.e)
                np.savetxt(f, np.array(out)[np.newaxis, :], fmt="%.8f")

            # Save phi, wavefunctions, and density
            io.save(phi_out_file.format(vset=vset, n=n), {"phi": d.phi})
            wf_dict = {f"WF_{i}": np.abs(set_device.eigenfunctions[:, i])**2 for i in range(set_device.eigenfunctions.shape[1])}
            io.save(wf_vtu_file.format(vset=vset, n=n), wf_dict, set_device.mesh)

            plot_slice(set_device.mesh, np.abs(set_device.eigenfunctions[:, 0])**2, 
                    normal=normal, origin=origin, 
                    cb_axis_label="Ground-state density [$e$ / m$^3$]",
                    path=wf_slice_file.format(vset=vset, n=n), show_figure=False)
            plot_slice(d.mesh, d.rho / ct.e, normal=normal, origin=origin, 
                    cb_axis_label="$\\rho$ [$e$ / m$^3$]", 
                    path=rho_slice_file.format(vset=vset, n=n), show_figure=False)


            # Compute many-body ground-state energy
            mb_E[j, i] = np.sum(g * set_device.population_factors * set_device.energies)
            print("  Many-body ground-state energy (eV): "f"{mb_E[j, i] / ct.e:.10f}")

            if i > 0:
                # Compute chemical potential
                Mu[j, i-1] = (mb_E[j, i] - mb_E[j, i-1])
                print(f"  Chemical potentials (eV): {Mu[j, i-1] / ct.e}")
                with open(mu_file, "a") as f:
                    np.savetxt(f, np.array([[vset, n, Mu[j, i-1] / ct.e]]), fmt="%.10f")

    return Mu