diamond internals

This file documents the internals of the diamond program. If you want to modify the program, this may be worth reading before the sources. I assume you have read the article Solving rate equations for electrons tunneling via discrete quantum states which explains what all this is about.

This program is intended to simulate the I-V curves measured on a single electron transistor in the regime where tunneling takes place through individual electron quantum states. It reads a configuration file containing the description of the device as well as the experimental conditions. The latter can be overriden by command-line options. The I-V curve is output to stdout as numbers in ASCII form.

Program structure

Files

The program consists of four C files and associated headers:

memory.[ch]: A few wrappers around malloc(), calloc() and free().
parse.[ch]: Code to parse the configuration file. Exports the parse() function.
diamond.[ch]: This is where the calculation is done. Exports the get_steady_current() function.
main.c: Parses the command line and wraps everything together: calls parse() to get the device description and makes a loop around get_steady_current() to output the current as a function of the parameter(s) chosen by the user.

The dependences are:

memory is used by diamond and main
parse is used by main
diamond is used by main

In the following, only the contents of diamond.[ch] will be discussed.

Inside diamond.c

The exported interface of diamond.c is mainly made by the following prototype:

double get_steady_current(device_t *device,
    double temperature, double e_r, double e_l);

where the parameters are respectively:

a description of the device,
the temperature in K,
the right lead Fermi energy,
the left lead Fermi energy,

and the return value is the steady-state current in pA. Here is a (slightly) simplified version of this function:

double get_steady_current(device_t *device,
    double temperature, double e_r, double e_l)
{
    problem_t *pb;          /* rate equation problem */
    double *population;     /* steady state population */
    double current;

    pb = generate_pb(device, temperature, e_r, e_l);
    population = get_steady_state(pb);
    current = get_current(population, pb);
    return current;
}

As can be seen above, the work is done in three steps:

generate_pb() constructs a rate equation problem from the input parameters. This problem is described by a structure of type problem_t.
get_steady_state() solves the rate equation problem for the steady-state population.
get_current() uses this solution, together with the original problem, to get the actual current.

Data representation

Device description

The device is described by the following structure:

typedef struct {
    int n_l;                  /* number of levels */
    double *level_pos;        /* level positions in meV */
    double *g_r, *g_l;        /* right and left gammas in MHz */
    int n_e;                  /* number of electrons */
    double **interactions;    /* interactions matrix in meV */
    double interactions_max;  /* maximum dispacement of a level */
    double r_b, r_g;          /* bias and gate prefactors (unused) */
} device_t;

which is filled from the configuration file. level_pos, g_r and g_l are arrays of length n_l. interactions is an n_l * n_l matrix. interactions_max is an estimate of the maximum displacement a level can experience. It is used to discard levels that are too far above or below the Fermi energy of both leads. r_b and r_g are numerical prefactors used in main.c to get the Fermi energies from the bias and gate voltages. They are ignored in diamond.c.

Many-body states

The problem_t structure has basically to describe the set of available many-body states and the transition matrix between these states. In order to get rid of the factors 2 that come from spin degeneracy, we will consider single-electron states wich differ only by the spin of the electrons as different states. Then, there are twice as many single-electron states as energy leves. A many-body state is described by a sequence of bits, each giving the occupancy of a single-electron state. This is just a binary number, as in:

typedef unsigned long config_t;

Here config_t stands for ``electronic configuration''. From here on, we will call configurations the many-body states. The single-electron states will just be called states. The type definition above sets the maximum number of levels to 16 assuming 32 bits longs.

The set of availables configurations is represented by the following fields of the problem_t structure:

    int n_u, n_c;           /* nb of uncharged and charged configs */
    config_t *conf;         /* available electronic configurations */
    int *reverse;           /* reversed conf[] */

Uncharged configurations are those with n_e electrons, while charged ones have n_e+1 electrons. conf is an array of length n_u+n_c whith the n_u first elements being uncharged and the n_c next being charged. The array reverse is used tho find the position in conf of a given configuration. In other words, for any given element conf[i] of conf, we have reverse[conf[i]]=i. It should be noticed that reverse is an array of length 2ⁿ⁺¹ when there are n levels available.

Rate equation problem

It is impractical to represent the transtion matrix as a plain matrix because of its potential size. For example, there are close to 3.5*10⁴ configurations for 10 or 11 electrons in 10 levels. A square matrix of this size would require close to 1 TB of memory. Fortunately, most of its coefficients would be zeros.

Since we are neglecting internal relaxation, cotunneling and other high order processes, the only non-zero coefficients in the transition matrix are the diagonal coefficients and those coupling configurations that differ by only one electron. We will call neighbors such configurations. All uncharged configurations have the same number of neighbors n_n, which are all charged configurations. This number is just the number of holes in an uncharged configuration. We use a 2D array neigh_i to give the index in conf of the neighbors of the uncharged configurations: the k^th neighbor of configuration conf[i] is conf[neigh_i[i][k]]. Then these indices i and k are used to identify all the possible transitions.

Below is the structure describing the rate equation problem:

typedef struct {
    int l_l, l_h;           /* lower and higher relevant levels */
    double **g_i, **g_o;    /* tunneling-in and tunneling-out rates */
    double **g_f, **g_b;    /* forward and backward rates in right lead */
    config_t *conf;         /* available electronic configurations */
    int *reverse;           /* reversed conf[] */
    int n_e;                /* number of electrons */
    int n_u, n_c;           /* nb of uncharged and charged configs */
    int n_n;                /* nb of neighbors for each uncharged state */
    int **neigh_i;          /* index of neighbor of uncharged state */
    int **neigh_l;          /* level changed in neighbor, XXX: remove. */
} problem_t;

In this structure the four g_* fields are n_u * n_n matrices, where g_*[i][k] represents some transition rate between configurations conf[i] and conf[neigh_i[i][k]]. These fields represent the following:

g_i: tunneling-in rate, i.e. g_i[i][k] is the rate of the conf[i] -> conf[neigh_i[i][k]] transition
g_o: tunneling-out rate, i.e. g_o[i][k] is the rate of the conf[neigh_i[i][k]] -> conf[i] transition
g_f: forward current in the right lead, i.e. the contribution of the right lead to the g_i[i][k] transition rate.
g_b: backward current in the right lead, i.e. the contribution of the right lead to the g_o[i][k] transition rate.

This structure allows to write the product of the transition matrix with a probability vector p in a fairly efficient way. This is done by theget_rates() function. Here is its inner loop (pb is the problem descriptoin and j an int, v will hold the product):

            /* Get index of neighbor. */
            j = pb->neigh_i[i][k];

            /* Charging transition. */
            v[j] += pb->g_i[i][k] * p[i];
            v[i] -= pb->g_i[i][k] * p[i];

            /* Uncharging transisition. */
            v[i] += pb->g_o[i][k] * p[j];
            v[j] -= pb->g_o[i][k] * p[j];

The current is calculated by get_current() in a very similar way using g_f and g_b.

Last updated: 2001-07-20

Edgar Bonet
edgar@edgar-bonet.org.