Speeding up device simulations

If you run a large number of simulations, the speed at which each one completes matters. This page explains how to make OghmaNano's electrical solver run faster. Because the useful settings only make sense once you understand what the solver is actually doing, the first half of the page covers the foundations of the model — the equations that are solved and how they are solved — and the second half shows how to tune the solver to shorten run times. The physics is unchanged by any of these settings; they affect only how quickly the numerical solution is reached.

1. The foundations of the model

Before looking at how to speed a simulation up, it helps to know what OghmaNano is computing. To account for charge transport in a device, OghmaNano solves the basic drift–diffusion equations together with Poisson's equation (Gauss's law). These are solved on a finite-difference mesh across the device, so that charge movement from one contact to the other can be tracked in both the conduction (LUMO) and valence (HOMO) levels.

1.1 Charge transport across the device

Poisson's equation relates the electrostatic potential to the free charge, while the two current equations describe the electron and hole driving terms (drift plus diffusion):

These current densities are then substituted into the electron and hole continuity equations, which balance the flow of current against recombination, trapping and the time rate of change of the free-carrier populations:

Here \( \phi \) is the electrostatic potential, \( n \) and \( p \) the free electron and hole densities, \( \mu_e \) and \( \mu_h \) the mobilities, \( D_n \) and \( D_p \) the diffusion coefficients, \( E_c \) and \( E_v \) the band edges, \( R \) the recombination rate, and \( T \) the net rate at which carriers are captured into traps. The full derivation and sign conventions follow MacKenzie et al., J. Phys. Chem. C (doi:10.1021/jp200234m).

1.2 Charge trapping and recombination

Charge trapping and recombination are important in organic devices and in many other material systems, so they must be modelled explicitly. It is common to reduce Shockley–Read–Hall (SRH) recombination to a single steady-state expression, which is easy to solve. OghmaNano does not do this: because the goal is to solve the equations in the time domain, the full set of SRH capture/escape rate equations is solved at each mesh point on the finite-difference grid.

At each point the situation shown in ?? is solved. The parabolas at the top and bottom are the free-carrier bands at the electron and hole mobility edges. Descending from the electron trap edge \( E_e^d \) is an exponential distribution of electron traps (red); rising from the hole trap edge \( E_h^d \) is an exponential distribution of hole traps (blue). Every one of the discrete levels within these distributions is an individual trap, and for each one the full capture-and-escape rate equations are solved. The change in the trapped-electron population of a level is the balance of four rates:

where \( r_1 \) is electron capture, \( r_2 \) electron emission, \( r_3 \) hole capture and \( r_4 \) hole emission, each following the usual capture/emission form set by the trap density, the carrier density, the capture cross-section and the thermal velocity (see the reference above). Detailed balance is maintained across all levels. Summing the net capture over the trap distribution gives the free-carrier recombination rate:

2. How the equations are solved: Newton's method

All of the equations above are solved with Newton's method, published by Isaac Newton in 1671 in Method of Fluxions. It is one of the most efficient and robust methods available for this type of coupled equation set: if an equation cannot be solved by hand, Newton's method will almost always find an answer.

2.1 Roots and the error function

Start with something easy. For what values of x does the following equation equal zero?

By inspection the answer is \( x = 3 \) and \( x = 40 \). These are the roots of the equation. The same problem can be written as a function,

and if \( f(x) \) is plotted against \( x \), the roots are simply the places where the curve crosses the axis, i.e. where \( f(x) = 0 \). Thinking of the equation as an error function that is driven to zero is the key idea. In the device model, \( x \) is replaced by the physical unknowns — the potential \( \phi \), and the electron and hole quasi-Fermi levels \( E_{Fn} \) and \( E_{Fp} \) — and the values of those variables are adjusted until \( f(x) = 0 \).

Roots are only obvious by eye for simple equations. Consider instead

Plotting this to find its roots would be a lot of effort, and the device equations are far worse. This is exactly where Newton's method earns its place.

2.2 The method, step by step

Take a concrete example, \( f(x) = x^2 + 10x - 3 \). Newton's method proceeds as follows:

1. Differentiate the equation with respect to the variable being solved for: \( f'(x) = 2x + 10 \).
2. Guess the root. Any value will do — take \( x = 10 \). It is a poor guess (\( f(10) = 197 \), nowhere near zero), but that does not matter.
3. Evaluate the function and its derivative at the guess: \( f(10) = 197 \) and \( f'(10) = 30 \).
4. Compute a better guess using Newton's update:

Substituting the values gives a better estimate, \( x_{+1} = 10 - 197/30 = 3.433 \). The procedure is then simply repeated: feed the improved guess back in, evaluate again (\( f(3.433) = 43.115 \), \( f'(3.433) = 16.866 \)), and obtain a still-better estimate of \( 0.877 \), and so on. Each pass drives the error function closer to zero and the estimate closer to the true root.

2.3 A worked example

The whole method is only a few lines of code. The loop below implements it for \( f(x) = x^2 + 10x - 3 \), printing the current estimate and the error at each pass:

x = 10;                       % our guess for the answer
for i = 1:10                  % loop ten times
    y  = x*x + 10*x - 3;      % f(x)
    dy = 2*x + 10;            % f'(x)
    new_x = x - (y/dy);       % calculate a better x
    x = new_x;                % update x
    sprintf("%f %f", x, y)    % print the estimate and the error
end                           % repeat

The output is instructive. The estimate moves \( 10 \rightarrow 3.433 \rightarrow 0.877 \rightarrow 0.321 \rightarrow 0.292 \rightarrow 0.291503 \), and the error drops \( 197 \rightarrow 43 \rightarrow 6.5 \rightarrow 0.31 \rightarrow 0.00085 \rightarrow 0 \). Notice how the number of correct decimal places roughly doubles each pass. This is called super-linear convergence, and it is what makes Newton's method such an efficient way to find the root of an equation.

2.4 The same idea in matrix form

The device is not a single equation but a large coupled set — the potential, the electron and hole Fermi levels and the trap occupations at every mesh point. Written in matrix form, Newton's method becomes:

\( [J] \) is the Jacobian, the matrix of derivatives of all the device equations (Poisson, drift–diffusion, SRH) with respect to all the unknowns. \( [f(x)] \) is the error in those equations for the current guess, and \( [\Delta x] \) is the resulting update — the improvement to the field, the Fermi levels and the trap populations across the whole device. The update is added to the current solution,

and the process is iterated until \( f(x) \approx 0 \), at which point the device equations are solved. Two points carry directly into the tuning that follows: \( f(x) \) is the quantity we are driving to zero, so it is the natural measure of convergence; and the better the initial guess, the closer the solver starts to the answer, and the fewer iterations it needs.

3. Speeding up a simulation in practice

With the foundations in place, the practical settings make sense. The example used here is a two-dimensional field-effect transistor (a 2D device, with all traps included), but the same ideas apply to any simulation.

3.1 Turn off writing to disk

The single easiest change is to stop writing so much data to disk. Writing files is slow and, for a simulation that produces a great deal of output, it can dominate the run time. OghmaNano provides a control that minimises the amount of data written to disk in one click; the individual output-file options (for example the dumping of dynamic time-domain data) can also be toggled on and off individually in the output settings, so you can keep only the files you actually need. For the fastest possible run, write as little to disk as you can.

3.2 Watching the solver converge

To tune the solver you first need to see it working. OghmaNano can print the Newton (current) error for each iteration and write the convergence history to a file, rather than printing a single line per voltage step. Enable the Newton-error output and the convergence log, then re-run.

The first simulation step is special. When the model starts it has no good guess for the device state — it does not yet know the potential or the Fermi levels — so it begins a long way from the solution. In ?? the error starts enormous (around \( 10^{15} \)) and is driven down over many iterations until it reaches a noisy floor near \( 10^{-7} \) to \( 10^{-9} \), at which point the equations are essentially solved. Because this first solve is so far from the answer, two things are done deliberately: the Newton step is heavily clamped (see below) so the solver converges slowly and safely rather than jumping somewhere that might crash, and it is allowed a large number of iterations — typically around a thousand — to give it a proper chance to settle.

After the equilibrium solve, the model steps through the requested voltages. Each voltage step moves the device slightly away from the previous solution, so the error jumps up and is then very quickly driven back down — the sequence of spikes on the right of the plot. Because each step starts close to the last answer, these solves converge rapidly.

3.3 Tuning the solver

The convergence plot is exactly the error function \( f(x) \) from Section 2, and the goal of tuning is to make it shorter — to reach convergence in fewer iterations. The main controls in the solver configuration are:

• Maximum iterations on the first step. If the first solve has clearly reached the noise floor after, say, 400 iterations, there is no value in running the full thousand. Reducing this limit (for example from 1000 to 600) simply trims the flat tail where the solver is already converged. Each removed iteration is real computational time saved.
• Target error on the first step. The first solve is often set to reach a very small error such as \( 10^{-9} \). If the result is already well converged at \( 10^{-7} \), relax the target: stopping at \( 10^{-7} \) instead of chasing \( 10^{-9} \) saves the last, slowest part of the descent.
• Clamping. The clamp is the fraction of the full Newton step \( [\Delta x] \) the solver is allowed to take. A full step divides straight through by the Jacobian; clamping deliberately takes only part of it — if Newton asks to move a variable by 0.5, a clamp might allow only 0.1. Small clamps converge more slowly but are far less likely to send the solver somewhere it crashes; larger clamps converge faster but more riskily. The first step usually uses a small, conservative clamp (around 0.1); raising it (0.2, 0.3, and higher on well-behaved structures) sharpens the initial descent noticeably. The clamp for the remaining steps can likewise be increased from its default, though pushing it too far risks non-convergence.
• Minimum number of Newton iterations. This forces the solver to run at least a set number of times before it is allowed to exit, which helps it climb out of any false minimum it might otherwise stop in. A small value (a few iterations) is useful insurance; setting it too high just wastes time bumping along the noise floor after the solution is already found.
• Minimum solver error and maximum iterations (all other steps). These set the target accuracy and iteration ceiling for the per-voltage solves. They behave just like their first-step counterparts and can be tightened once you know how quickly your structure converges.