Simulations of Schrödinger-Poisson converge self-consistently, and almost automatically, thanks to an algorithm proposed some years ago by Alex Trellakis, one of our talented developers. This algorithm was implemented in nextnano++ and have been very successful for several devices.

However, when the current equation is included in this system, the convergence to the solution becomes a challenge, due to the nature of this equation. For some devices, the system of equations becomes very unstable and a certain ability to reach the convergence is required.

Especially for systems where the carrier density fluctuates from large values to almost zero in certain regions or interfaces, the process of obtaining convergence becomes more critical and acting in a strategic way is very helpful.

Setting the input file for performing self-consistent current-Schrödinger-Poisson computations

Self-consistent current-Schrödinger-Poisson computations can be specified in the section run{ } of the input file, through the statements

  • current_poisson{}

  • quantum_current_poisson{}

The first statement is mandatory, and it provides a first estimate of the electrostatic potential and the (quasi-)Fermi levels, even before including the quantum calculations to the system.

In principle, this is the minimum information required to start the simulations. All numerical parameters are adjusted automatically internally in the code until the solution is found or the maximum number of iterations is reached. Unfortunately, given the huge variety of devices the program can simulate, universal parameters are not possible to be predicted in advance. For this reason, in order to give to the user more control of the convergence process, some parameters can optionally be specified within the subsection quantum_current_poisson{}. Some examples are the following parameters: alpha_Fermi, residual, residual_fermi and iterations.

It is not our purpose to describe each of these parameters in this document, but to provide some guidance how to control the numerical process with the minimum effort as possible. The list of all parameters, its description, range of values and default can be found on the section run.

Talking about convergence

Before proceeding it is important to discuss what the expression “to get convergence” means. Actually, nextnano++ has to solve six groups of equations: current (also called, continuity) equation, Schrödinger equation and Poisson equation for electrons and for holes. As default the values of carrier densities, Fermi levels and potential are kept iteratively consistent from one step to the other. Internally the program computes for each equation a so-called cost function, that represents a metric of how close the obtained solutions are close to the “exact” one. For example, one way the cost function can be defined is by the difference of left and the right of each equation. Then, after each iteration the results of the cost function are called residuals.

Getting convergence means to find the conditions that minimize the cost functions. A good analogy of this process is the task of finding the deepest location of a valley in a mountain chain. In order to reach this valley, having some strategy concerning the necessary moves in some direction can reduce the time and the number of steps to conclude this task. If each step is too large, we can overfly the valley, if it is too small, we can take a long time to reach it. This is the role of the alpha_Fermi parameter in the current_poisson and quantum_current_poisson solvers: large values of alpha_Fermi can make the minimum invisible, and if it is too small can take a long time for simulations. Additionally, especially when the value of alpha_Fermi is small, it is possible that the number of iterations, given by the parameter iterations, is not enough to reach this minimum.

This analogy with a mountain chain is actually very simplistic, because the program deals with finding a minimum of cost functions in a multi-dimensional space and a non-linear system of equations, which makes this task more complex and, for this reason, provides more accurate results than any analytical model.

Keeping this in mind, setting the right parameters is usually an iterative process. One procedure that can be used for reducing the simulation time is by displaying the results “on-the-fly” within our graphical interface (nextnanomat) in two different ways.

The first method is through the numerical values displayed in the “Simulation” tab of the graphical interface (nextnanomat). The evolution of the residuals is printed out as soon they are computed. If, after certain time, some of the residuals are not reducing from one iteration to the next one, it is recommended to stop the simulation and restart a new one with different parameters.

The second method is by plotting the files interation_current_poisson.dat and iteration_quantum_current_poisson.dat. By default, these files are generated automatically by the program, unless “output_log = no” is specified in current_poisson{} and quantum_current_poisson{} subsections. They can be displayed in the browser menu for the “Output” tab of nextnanomat. As in the previous method, if the residuals are reducing too slow, it is recommended to restart a new simulation that can accelerate the process.

Getting some intuition…

As mentioned above, depending on the nature of the device and the specific operation conditions (temperature or bias), it is necessary to guide the tool to get convergence. Let us see some practical examples.

Here we will illustrate how the evolution of the residuals in a current-Schrödinger-Poisson can evolve during the convergence process for two different devices. The images correspond to the plot of the data from interation_current_poisson.dat, and iteration_quantum_current_poisson.dat files, that can be found in the output folder of the simulation.

Figure 1 corresponds to the residual evolution of a system that converges faster: all residuals drop around one order of magnitude every ten iterations. The default parameters within the code brings the system almost automatically to the minimum of the residuals.

In contrast, Figure 2 shows the final result for a different device after the system gets convergence. In this case, in the input file were especified that residual_fermi is equal to \(10^{-7}\) eV and residual (density) as \(10^5/cm^3\). The value of alpha_Fermi is 0.01. Although it was specified a total of 2000 iterations, the convergence was achieved in around 400 steps. It is important to notice that only after 200 iterations the system starts reducing the residuals in several orders of magnitude.

For some devices, setting the values of alpha_iterations and alpha_scale can result in a better performance. The value of alpha_iterations is related to the moment where the alpha_Fermi shall start to gradually reduce, and the value alpha_scale is the rate of reduction between two successive iterations. There is no rule for the direction they should be changed. It is necessary to test some cases and look at the effect on the residuals.


Figure 1 Residual evolution for a system exhibiting quick convergence.


Figure 2 Residual evolution for a different system than the one from Figure 1. In the input file were specified residual_fermi = 10^-7 eV, residual (density) = 10^5 /cm^3, and iterations = 2000. The alpha_Fermi parameter was set to 0.01.

Sometimes the number of iterations is not enough to reach the convergence. Figure 3 corresponds to the plot of results for the same system of the Figure 2, simulated with only 150 iterations. As it was shown in the previous figure, only after 180 iterations the residuals start to decrease. In this kind of simulations, there are no criteria for knowing at which point this will happen: it requires experience or can be done by trial and error.


Figure 3 Residual evolution of the same system of Figure 2, exhibiting a pseudo-non-convergence behavior. Specifications in the input file: residual_fermi = 10^-7 eV, residual (density) of 10^5 /cm^3, and iterations = 150. The value of alpha_Fermi is 0.01.

A pseudo-non-convergence can also happen when small residuals are specified in the input file. Returning to the Figure 2 it can be observed that, choosing residual_fermi as 10-^10 eV would probably result in a non-convergence: the residual_fermi does not decrease at a high rate after 350 iterations. Then, increasing the number of iterations in this case would not solve the problem.

Another situation is when the value of alpha_Fermi is too small: it looks like the residuals do not reduce, like in Figure 4. In this example, alpha_Fermi was reduced from 0.01 (value used for figures 2 and 3) to 0.0001, and after 2000 iterations the system does not converge. Here we used the same system of the previous two images. There are other patterns, but here the most relevant one is shown.


Figure 4 Residual evolution for a system exhibiting pseudo-non-convergence. Specifications in the input file: residual_fermi =10^-7 eV, residual (density) = 10^5 /cm^3, and iterations = 2000. The alpha_Fermi parameter was reduced to 0.0001.

Sweeping parameters

It is very common to use a sweep of specific variables within the input file, for example bias or any other user defined parameter.

It is important to have in mind that any change in the input file is equivalent to a simulation of a new system (for example when modifying doping), or the operation condition (temperature or bias). There is no mathematical reason that the solutions of two systems shall be similar. In other words, it is not expected that all solutions using different conditions will converge under the same criteria, for the entire range of variation of the sweep parameters. Eventually, for example, a sweep of bias from 0 to 8 Volts can use the same parameters for the whole simulation, but this is not the most common case.

A good strategy is to start the sweep of the parameters and verify at which value the solution does not longer converge. For saving time it is recommended to split the range of variation in two parts, and to follow the simulation only using the values of the parameter (for example, bias) that have still not converged. Trying to make the solution converge for a wide range of values for the sweep variable, using with a unique set of residuals and alpha_fermi, can become a very hard task, without the recommended range splitting.

… and when nothing works

Our concern, in the development of our code, is to make it as accurate and fast as possible. Some simulations can be performed in a simple notebook, especially for 1D simulations.

Unfortunately, for some devices under specific conditions, making the system of Current-Schrödinger-Poisson converge in few iterations is a very specialized and time-consuming task. Observing the needs of our customers, nextnano is offering our customers the opportunity to perform this task on demand. Please, consult our schedules and fees when an extra assistance is required. Our experts in simulation can assist you to boost your project!

Copyright © 2021 by nextnano GmbH

This document is copywritten by nextnano and it is distributed for free to nextnano customers. The company does not allow sharing, publishing or distributing by any means without its written permission.

Edited by Maria Cecilia da Silva Figueira.