SEMITIP V6, Technical Manual

Technical instructions are provided below for SEMITIP version 6, in the sections: Additional information on the theory behind the SEMITIP package can be found in Refs. [1]-[6] and in the documents:


The SEMITIP version 6 software package contains a set of programs for computing the electrostatic potential and resulting tunnel current between a metallic tip and a semiconductor sample. An brief introduction to the package is provided at the main SEMITIP V6 Web Page, including a discussion of various dimensionalities, geometries, and types of tunnel current computation that can be handled. A different program is employed for each dimensionality, geometry, or type of current computation. Links to those programs, including documentation for each, are included on the main SEMITIP V6 Web Page. The purpose of this Technical Manual is to describe various methods employed in the algorithms, coding, or usage of the SEMITIP package that are common to a range of the programs.

As mentioned above the SEMITIP package permits computations for many different dimensionalities, geometries, and types of tunnel current computations. The goal of the package is to both provide users with the capability to easily perform computations using existing programs suitable for their problem, and to enable them to develop new routines or programs for handling situations not covered by presently existing programs. This first goal could in principle be handled by a single master program that handled every possible type of situation. That code would, however, be very unwieldy since by its nature it would accomodate many different situations. For development of new code by users, generally there is some specific situation (i.e. dimensionality, geometry, type of current computation, etc.) that is being considered. Then, it is very much easier for a user to modify existing code that deals with that specific situation rather than dealing with code that covers every possible situation. This is the reason that separate programs have been developed to deal with the various dimensionalities, geometries, and types of tunnel current computation that are covered within the SEMITIP package.

Input Parameters

All of the programs in the SEMITIP package take their input from a file named FORT.9 (i.e. with a different FORT.9 file for each program). Comments within each of those files specify the meaning of each parameter, with additional information for certain parameters provided below. The line numbers in the description below are applicable to the Uni2 program for parameters having to do with computation of potential, to the UniInt2 program for parameters having to do with computation of current by integration of the Schrödinger equation, to the UniIntSC2 program for parameters having to do with self-consistency when computing the current by integration of the Schrödinger equation (only relevant for situations of inversion or accumulation), and to the UniPlane3 program for parameters having to do with computing the current using a plane wave expansion. However, the relevant line numbers for the FORT.9 files of other programs can be easily deduced from those listed below since the parameters within every input file are clearly labeled within the file itself.

Parameters for computation of potential: (line numbers according to FORT.9 file of Uni2)

Lines 7-18 refer to parameters defining the semiconductor. For an inhomogeneous semiconductor containing various regions of different doping or various surface area with different surface charge densities, see section on Inhomogeneous Geometry. Lines 21-28 refer to parameters that define the distribution in energy of surface states on the semiconductor. By default two types of distributions can be input, and the charge densities from these are summed together in the program for computing the total charge density (physically, many surface might have intrinsic and extrinsic surface states, the former from the states of the intrinsic surface atoms and the latter from defects on the surface, and hence the use of two distributions). The distributions themselves are defined in a user-defined function SIG which is part of the main program (for further discussion see section on User-defined Functions). The parameters below corresponding to the default form of the SIG function. By default two types of distributions are possible: a uniform distribution of states and a distribution containing two Gaussians. In the first case the parameters are just the density of states and the charge neutrality level [i.e. that separates states of acceptor character (negatively when occupied by an electron and neutral when empty) from those with donor character (neutral when occupied by an electron and positive when empty)]. To use this type of uniform distribution, then the value of the FWHM parameter should be zero (and the centroid energy is not used by the program, although it still should be present in the FORT.9 parameter list). For the Gaussian type distribution, nonzero values for the FWHM and centroid energies are supplied. The former is the full-width-at-half-maximum for each Gaussian and the latter is the position of the peak of each Gaussion relative to the charge neutrality level (there's one Gaussion below that level and one above it). One additional issue regarding surface states is whether or not the distributions change (e.g. go to zero) for energies with the bulk valence or conduction band. This constraint can applied by a line of the code within the SIG routine, and for the default versions of SIG there is no change of the distributions when the energies fall within the bulk bands.

Lines 30-36 refer to parameters that determine details of the finite-difference solution for the potential. This solution is performed first on an a specified number of grid points. Then, the number of grid points is doubled (and their spacing halved), and the computation performed on that doubled grid (using as an initial condition the solution from the previous grid). This procedure of doubling the number of grid points in each coordinate direction is referred to as "scaling" of the solution, and the scaling procedure is repeated a specified number of times. The finite-difference procedure itself is an iterative one, and within each scaling step it is performed until the change in Pot0 (the potential at the point on the semiconductor surface opposite the tip apex) is less than a specified convergence parameter value for both the present iteration and the previous one, or the number of iterations exceeds a specified maximum number. Concerning the initial spacing between grid points, this is straightforward for the points in the vacuum, but in the radial direction or in the direction into the semiconductor some estimation of the initial spacing must be made. This estimation of made based on values of semiconductor doping and tip radius as follows: First a 1D depletion length of the semiconductor is computed, sqrt(2 epsilon ΔV / e2 N ) where epsilon is the dielectric constant and ΔV corresponds to the tip-sample voltage or to 1 V, whichever is greater. The grid size is initially assigned the value of the tip radius. If the radius of a protrusion on the end of the tip is nonzero, then the grid size is assigned to that value or to its previous value, whichever is less. Then, the grid size is assigned to the above 1D estimate of depletion length divided by the number of grid points, or to its previous value, whichever is less. Finally, this value for the initial grid size is multiplied by the factor on line 33 (e.g. to reduce its value further, and hence achieve a finer grid). The grid size thus determined is used for both the radial direction and the z-direction into the semiconductor, except if the multiplier parameter on line 33 has a value <= 0. In that case, the following two additional lines provide the values of the grid size in the radial direction and the z-direction (or, for a 1D computation, just one additional line is used for providing the grid spacing in the z-direction). The resulting grid sizes are used to compute the spatially varying spacing between grid points using the algorithm described in the section on Spatial Coordinates.

Parameters for computation of tunnel current by integration of Schrödinger's equation: (line numbers according to FORT.9 file of UniInt2)

Computation of Bulk Charge Densities

Bulk charge densities are computed in the effective-mass approximation, using equations as specified in Ref. 1. These charge densities are evaluated repeatedly within the programs - at each grid point, each iteration of the solution, and each 1D search step to solve the nonlinear aspect of Poisson's equation. To reduce the computation time needed for these evaluations, a table of charge density values to computed at the start of each program (or multiple tables are computed, for an inhomogeneous semiconductor). The maximum possible number of elements in this table is specified by NEDIM in the PARAMETER statement within each program, and the actual number of element used is specified in the FORT.9 input file (line 37). The program uses a particular algorithm to decide on the energy bounds of this this table, and if a charge densities value is needed that falls outside of these bounds then that charge density is evaluated from the defining equations.

Computation of Surface Charge Densities

The presence of electronic surface states give rise to surface charge densities. There are in general two broad categories of surface states that can exist on the semiconductor, one being intrinsic states from the dangling bonds that exist within each unit cell of the structure, and the second being extinsic states that are sparsely distributed over the surface and arise from defects such as adsorbates or surface steps. For surfaces such as GaAs(110) that do not have any intrinsic states within the band gap, tip-induced band bending plays a particularly large role. Such surface can still have extrinsic surface states, which will affect this band bending. (Also, even for GaAs(110), intrinsic surface states located within the conduction band play an important role for situations of electron accumulation, as discussed in Ref. [5]). By default, the SEMITIP programs allow for input of two type of surface charge densities. Both of these densities are treated identically within the program, with the reason for allowing two inputs being to treat situations such as GaAs(110) for which both intrinsic and extrinsic states can play a significant role in determining the band bending.

Surface charge densities are defined by user-specified distribution(s), defined in the routine SIG. In the default form of those routines, a distribution that is uniform in energy, or one consisting of two Gaussian functions, are provided, but any other function can also be specified. (In principle, a function that varies very quickly with energy may be difficult for the program to handle, but no such limitation has yet been encountered). In addition to the energy distribution, the charge neutrality level for the distribution must also be specified, which is the energy above which the states are negative when filled and neutral when empty (i.e. acceptor like, as for conduction band states) and below which they are neutral when filled and positive when empty (i.e. donor like, as for valence band states). In analogy to the procedure for bulk charge densities, a table of surface charge densities is constructed at the start of each program. For the surface charge densities in particular this table assumes a zero temperature form for their occupation. Computations using occupations computed at nonzero temperature are possible (using the switch on line 29 of the FORT.9 input file), but in that case the table is not used and the occupations are always computed from the defining equations, thus requiring considerably more computation time. It is rare that the explicit inclusion of the temperature dependence of the occupation significantly affects the results, i.e. to an extent beyond what could be accomplished by a slight shift in one of the surface state parameters (such as the charge neutrality level), so that in most cases this temperature dependence need not be included in the computation.

It is important to realize that SEMITIP, in its present form, does not allow for any computations of tunnel current associated with the defined distributions of surface states. Rather, the definied charge density of the surface states is used only in computing the electrostatic potential, whereas the tunnel current is based entirely on bulk bands (possibly including localized regions of the potential as occur e.g. for a buried quantum dot). In some future version of SEMITIP it is possible that currents due to surface states could be included; the surface states in that case would be defined by a surface band structure, and the charge density for that band of states would be computed by an integral of the band structure (rather than be a user-defined density-of-states, as presently used).


A computation of electrostatic potential is said to be self-consistent if the solution for the potential contains full and proper dependence on the electronic states of the problem, with the electronic states themselves depending on the potential. The default mode of charge density computation within SEMITIP is a semi-classical one, in which defined bands of bulk or surface states are shifted in accordance with the electrostatic potential and they are occuppied according to a fixed (known) Fermi energy. For a semiconductor in depletion, this type of computation is automatically (trivially) self-consistent, since the only relevant charge density of the problem is that due to the ionized dopants and this charge density is accurately treated in the semi-classical approximation. Incorporating user-defined distributions of surface states into the solution for the potential is also self-consistent in the same trivial sense (assuming that those surface states are not modified by the application of the potential).

Nontrivial situations of self-consistency occur for accumulation or inversion in the semiconductor. These situations are not properly handled by a semi-classical treatment (at least not when only one or a few quantum states are occupied), but the UniIntSC1 or UniIntSC2 programs are specially designed to treat these cases. The quantum states are computed by the intcurr routine, i.e. the same routine that is used for computing tunneling currents, but now used for constructing charge densities are needed for the self-consistent computation. The computation of the 3D charge densities is performed as a series of 1D computations at different r values, something that is appropriate only for a potential that varies slowly in the radial direction.

For self-consistent computations with an even more localized potential, e.g. as might occur at a quantum dot or near a point charge on a surface, the computation of the charge densities would have to employ the plane wave expansion method, as in UniPlane3. The result would be limited to a limited region about the point on the surface opposite the tip apex, i.e. a periodic repetition of this region. In any case, a program for handling such computations has not been developed to date.

Specification of Fixed Charge on Surface or in Bulk

A problem of considerable interest is the computation of the potential (and possibly the current) in the presence of a fixed charge, e.g. a point charge, on or near a surface. SEMITIP was not originally designed to handle fixed charges; rather, the parameters within the FORT.9 file refer everywhere to charge densities that vary in accordance with the Fermi energy. Fixed charges, whose magnitude does not vary with the Fermi energy, can still be handled in the program, but they are handled in a slightly "kludgy" manner making user-defined modifications to the RHOSURF or RHOBULK routines. See example 6 of Uni2 or example 2 of Uni3. In these examples, the location and magnitude of the fixed charge is explicitly defined with the RHOSURF or RHOBULK routine. An alternative method, of course, is for the user to modify the input stream of the FORT.9 file to include such parameters, and they then would be passed to the RHOBULK or RHOSURF routines through a user-defined COMMON block.

It is important to note that a number of limitations exist for computations involving fixed charges:

  1. First, the reader is reminded that the default manner in which SEMITIP handles all charge densities is a semi-classical one. This method is only valid for a slowly varying potential. Some situations that require explicit quantum computations are also handled within SEMITIP, such as semiconductor accumulation or inversion, but the package does not automatically handle all such situations. For a point charge on or near a surface, the electrostatic potential when the semiconductor is in depletion can be handled by the program (subject to some additional limitations described below), but situations of accumulation around the point charge (that is, occupation of localized states arising from the point charge as well as screening by extended states around the point charge) are not properly handled by the present package since it does not handle self-consistent computations with a plane wave basis. (Even if it did handle such computations, the limited spatial extent of any plane wave type computation would also be a significant restriction on the results).
  2. A second limitation in defining fixed charge on the surface or in the bulk has to do with the grid used in the finite-difference computation. As illustrated in example 6 of Uni2 and example 2 of Uni3, the spatial extent of the charge is defined by statements within the RHOSURF or RHOBULK routines. However, if that spatial extent does not precisely overlap with the areas or volumes subtended by the grid points of the computation, then the actual charge that is computed will differ from the intended charge. Hence, a sufficiently fine grid and/or an extended charge must be used so that this discrepency is not too large. For a charge located some distance from the tip position (central axis), it's no problem to model that charge as being extended over some region (since that effect of the charge is essentially the same as if it were localized at a single point). But, as the charge approaches the central axis, then this area or region over which it is definied must become smaller and smaller, with the grid size for the computation necessarily being correspondingly small. The user must ensure that these sizes for the spatial extent of a fixed charge and for the grid are appropriate to the particular situation being considered.

Solution of Poisson's Equation

Poisson's equation, including the boundary condition at the semiconductor surface and possible nonlinearity arising from bulk or surface charge densities, is solved according to the method described in Ref. [1] (also see footnote 36 of Ref. [2] and point 1 of the Internal Parameters section below).

Spatial Coordinates and Grids

The solution for the potential is accomplished using cylindrical coordinate in the semiconductor and generalized prolate spheroidal coordinates in the vacuum. The latter are chosen to precisely match the shape of the probe tip (not including any protrusion on the end of the tip), as described in Ref. [4] and in Coordinate System with accompanying Diagram. The spacing of the grid points in these coordinates is not uniform, but it increases as the radial distance r away from the central axis increases or as the distance z into the semiconductor increases. For computation of the tunnel current, this type of grid with quite large step sizes at large r or z is not suitable, and hence new grids with more nearly uniform spacing are formed. This process is done by the routine POTEXPAND for a tunnel current computation by integration of the Schrödinger equation or by POTPERIOD3 for a computation using the plane wave expansion method).

Computation of Tunnel Current

SEMITIP presently does not allow computation of tunnel currents from surface states. For the case of bulk states, two methods are available for computation of the tunnel current. In both, the Bardeen method together with the Tersoff-Hamann approximation for a sharp tip is used to write the current as a summation of the quantum states of the semiconductor (and the metallic probe tip). To obtain those states, the first method employs a "central axis approximation" in which numerical integration of the Schrödinger equation is performed only along the central axis of the problem. This method would be an exact solution for a planar geometry; it is an approximation for a nonplanar geometry but it still works quite well even for probe tips as sharp as 1 nm radius-of-curvature (see Ref. [6]). However, if regions of localized potential exist in the semiconductor, such as those due to a quantum dot, then that method in not applicable. Hence, the second available method is an expansion using plane waves (matched to decaying exponential in the vacuum), from which the solution of the Schrödinger equation is found by the usual eigenvalue method. Due to the significant computational demands of this plane wave method, it can be applied only over a limited region of space centered around the point on the semiconductor surface opposite the tip apex.

Inhomogeneous Geometry

A new feature in version 6 of SEMITIP is the ability to handle inhomogeneous semiconductors, with the inhomogeneity existing in the bulk or on the surface. Actually, even in the prior versions of SEMITIP it was possible to handle certain special types of inhomogeneity (such as the presence of fixed charges) by explicit modification of the RHOBULK or RHOSURF routines, as described above under Specification of Fixed Charge on Surface or in Bulk. But within version 6 a general method for handling inhomogeneity in many of the parameters describing the bulk or surface has been implemented. Such inhomogeneity does not represent any problem within a finite-difference type of computation that SEMITIP employs for the electrostatic potential, although whether or not a computation of tunnel current can truly handle an inhomogeneity depends very much on the range of the varying potential, as further discussed under Computation of Tunnel Current.

Inhomogeneities within the bulk of the semiconductor are referred to as different "regions" of the material, and inhomogeneities on the surface are referred to as different "areas". The user defines the spatial boundaries of the different regions or areas by using user-defined functions RHOBULK or RHOSURF, respectively, as described in the following section. Parameters specifying the properties of the different regions or areas (e.g. their doping or their surface state density) can, in most cases, be input in the normal fashion from the FORT.9 input file. The Mult1, Mult2, and Mult3 programs give examples of these procedures. Again, for specification of fixed charges within the bulk or on the surface the situation is different, as described under Specification of Fixed Charge on Surface or in Bulk.

In the present form of SEMITIP V6, computations of tunnel current for an inhomogeneous geometry are made using only a constant effective mass in the semiconductor. That mass is the one from the first defined region of the semiconductor. Thus, even if different masses for the various regions are defined within the FORT.9 input file, only the effective mass from region 1 is used in the computation of the current. (This limitation could be removed with some modest extension of the intcurr routine, but that extension has not yet been implemented).

When multiple regions are defined in the bulk material, then the locations of their respective valence band maxima can differ. Many of the quantities intput or output from SEMITIP (e.g. Fermi energy, charge neutrality level, etc.) use the valence band maximum as a zero of energy. In the situation with multiple regions, it is the valence band maximum of the region located at the origin of the coordinate system that is used as a zero of energy for all quantities, i.e. even those quantities associated with other regions.

User-defined Functions

A variety of user-defined functions are included in the same file as the main calling program. These functions allow definition of functions which cannot, in general, be specified simply by numbers in the FORT.9 input file. For example, a protrusion of arbitrary shape can be appended onto the end of the hyperbolic tip. Or, surface states of any arbitrary distribution in energy can be defined and their effect included in the computations of the potential. All user-defined functions are already included in the general purpose programs of the SEMITIP package, with default values for the functions as specified below:


The vast majority of SEMITIP is written using standard FORTRAN features, and it thus should be compatible with any FORTRAN platform. The programs are however designed for use under Gnu FORTRAN, and possibly a few compiler-specific FORTRAN features may be employed. For example:

COMMON blocks

COMMON blocks serve an important role in the SEMITIP programs by enabling the passing of parameters between the main programs and the user-defined functions, and in some cases to various supporting routines as well. The declarations and role of the various COMMON blocks is specified below. In situations when the block is only used between the main program and a user-defined function (i.e. it does not appear in any supporting routine), then modification by the user to accomodate some new situation is possible.

Internal Parameters

Nearly all the parameters that control the functioning of the SEMITIP programs can be accessed through either the FORT.9 input file, the user-defined functions, or the PARAMETER and COMMON statements described above. There are however a few additional parameters that are internal to certain routines, modification of which require editing the source code and then re-compiling and linking. These internal parameters include:
  1. A one-dimensional search is used as part of the finite-difference method for solving the nonlinear Poisson equation, as described in Ref. [1]. This search is performed to an accuracy in terms of energy of max(PotTip,1)/106, where max(PotTip,1) is the potential on the probe tip or 1 eV, whichever is greater. This precision parameter occurs several times within the semitip1, semitip2, and semitip3 routines. It is difficult to imagine a situation where this parameter would need to be changed, since a worse precision would not save a significant amount of computation time and it's unlikely that a higher precision would ever be needed (except for situations with very low temperature and a small energy scale, in which case possibly the precision of the entire computation would also require modification).
  2. The IBC parameter specifying the boundary condition for the potential computation is set to 0 (Dirichelet boundary condition) at the top of the semitip1, semitip2, and semitip3 routines. See Boundary Conditions for more discussion of this parameter. As discussed there, a value of 1 for Von Neumann boundary conditions can be used, but it's useful only for particular cases and the value of 0 is most suitable for general purposes.
  3. For both bulk and surface charge densities, tables of those values are constructed in the main program as discussed in the sections on Computation of Bulk Charge Densities and Computation of Surface Charge Densities. The energy bounds of those tables is determined by a particular algorithm within the main program. For some special situations charge densities might be required that are outside this range, although that's not a problem since in those cases the densities are evaluated from the defining equations. Nevertheless, if it is desired for some reason to change the bounds of the charge density tables, than that must be done by modifying the source code in the main program.
  4. For self-consistent computations, there's a parameter called INITCD within the UniIntSC1 and UniIntSC2 programs that controls the initialization of the array used for constructing the charge densities. By default this parameter to set to 1, which corresponds to the array being initialized to semi-classical values. An alternative is INITCD=0 which initializes the array to zero, but that hasn't proved useful in any computations to date.
  5. For computation of tunnel current, a two-band model for the band structure is not used. There are situation however where use of this model is appropriate, as discussed in Kane's Two-band Model. That document also explains how to modify the intcurr routine to enable a two-band type computation.

Changes relative to SEMITIP versions 1-5

Due to the considerable expansion in the capability of SEMITIP version 6 compared to versions 1-5, a number of changes to the naming conventions for programs and variables were needed. These changes included:

Additionally, in version 6, considerable repackaging of the routines occurred so as to enable the common usage of particular routines across many of the specific calling programs. The new program structure is described in the sections above.

For users upgrading from versions 4 or 5 to version 6, a few small changes in the input files (FORT.9) should be noted. In version 6, all input related to the computation and plotting of the potential is given before that for the computation of current. Hence, e.g. whereas in versions 4 or 5, specification of parameters for the contour plots occurred at the very end of the input file, those entries now occur at the end of the entries for the potential computation, and they are followed by the input parameters (if any) for the computation of the tunnel current.

One technical change occurred in version 6, namely, the introduction of a new method for computing derivatives for the variable-spacing grids that are employed in the semitip.f finite-difference routine. This method is further described in Computation of Derivatives for Variable Grid Spacing.


1. R. M. Feenstra, Electrostatic Potential for a Hyperbolic Probe Tip near a Semiconductor, J. Vac. Sci. Technol. B 21, 2080 (2003). For preprint, see
2. R. M. Feenstra, S. Gaan, G. Meyer, and K. H. Rieder, Low-temperature tunneling spectroscopy of Ge(111)c(2x8) surfaces, Phys. Rev. B 71, 125316 (2005). For preprint, see
3. R. M. Feenstra, Y. Dong, M. P. Semtsiv, and W. T. Masselink, Influence of Tip-induced Band Bending on Tunneling Spectra of Semiconductor Surfaces, Nanotechnology 18, 044015 (2007). For preprint, see
4. Y. Dong, R. M. Feenstra, M. P. Semtsiv and W. T. Masselink, Band Offsets of InGaP/GaAs Heterojunctions by Scanning Tunneling Spectroscopy, J. Appl. Phys. 103, 073704 (2008). For preprint, see
5. N. Ishida, K. Sueoka, and R. M. Feenstra, Influence of surface states on tunneling spectra of n-type GaAs(110) surfaces, Phys. Rev. B 80 075320 (2009). For reprint, see
6. S. Gaan, G. He, R. M. Feenstra, J. Walker, and E. Towe, Size, shape, composition, and electronic properties of InAs/GaAs quantum dots by scanning tunneling microscopy and spectroscopy, J. Appl. Phys. 108, 114315 (2010). For preprint, see