Refactoring legacy physics code with Boost

I recently came across an old undergraduate project of mine from the days when I thought I would do semiconductor physics as a career. It was a short C program for Monte Carlo Simulation of a particular electron scattering mechanism in GaAs (which should help you identify the timeframe). It was an academic exercise and didn’t quite work right, but I got an OK grade. Since that time, both software techniques and the skill level of the author have advanced, so I decided to give it another look.

Acoustic Phonon Scattering

Electrical resistance in semiconductor crystals comes from various kinds of imperfections or interactions that prevent carriers from taking a direct path. An impurity ion, for example, distorts the crystal lattice and electrons that encounter it can be deflected. These encounters are called scattering events.

The type of scattering addressed by the assignment stems from the natural vibrations in the crystal at temperatures above absolute zero. These cause localized expansions and contractions of the lattice that travel through the crystal. For modeling purposes these can be considered a type of particle: a phonon. These “particles” move through the lattice at the speed of sound and can collide with carriers, scattering them into new trajectories.

Phonon Scattering

Here we see one possible path of an electron, experiencing four phonon collisions as it moves through the crystal from left to right. In between each collision it accelerates (or decelerates) smoothly under the influence of the applied electric field $\mathbf{E}$ .

The accumulation of these interruptions results in an average drift velocity from which we calculate the critical parameter mobility:

$\begin{align*} \mu = \frac{|\boldsymbol{V_d}|}{|\boldsymbol{E}|} \end{align*}$

where $\boldsymbol{E}$ is the applied electric field and $\boldsymbol{V_d}$ is the measured average velocity.

Monte Carlo Simulation

Producing an estimate of mobility consists basically of running a large number of random simulated collisions and calculating the average velocity achieved by the electron, then dividing by the applied electric field. The steps for each event are described in more detail elsewhere but for our purposes they are:

Randomly select the time of the next scattering event
Calculate the average velocity of the period
Randomly select the new trajectory

After a large enough number of simulated collisions we should converge on an average velocity.

Boost.Math Constants

My first step after converting to C++ (a matter of a few minor fixes) was to replace hard-coded math constants like $\pi$ and $e$ with values from the Boost.Math numeric constants library. Although it’s true that most standard libraries supply M_PI and M_E, these are macros, with all their attendant problems. Instead I can use typed constants, like:

constexpr double two_pi   = boost::math::constants::two_pi<double>();
constexpr double     pi   = boost::math::constants::pi<double>();

Boost.Units Physical Constants

In addition to strongly typed physical quantities, the Boost.Units library supplies a rich set of commonly used physical quantities that I could use instead of floating point literals. Dividing by their units gives dimensionless numbers that I can directly substitute in my equations. For example:

using namespace boost::units;
constexpr double dirac   = si::constants::codata::hbar /
                           (si::joule * si::second);
constexpr double m_e     = si::constants::codata::m_e /
                           si::kilogram ;

gives me the Dirac constant $\hbar$ and the electron mass $m_{e}$ , as doubles.

Now I can rewrite my original equation relating the crystal momentum $\boldsymbol{k}$ to the velocity $\boldsymbol{v}$ :

$v_{x} = \underbrace{\frac{E_{x}\hbar}{m^{*}m_{0}}}_\mathtt{\text{vel_const}} k_{x}$

from

double const vel_const = (10000.0)*(6.63E-34)/((2.0*3.14159)*(0.063)*(9.11E-30));

constexpr double Efield = 10000.0;   // Applied field (V/cm)
constexpr double meff = 0.063; // effective electron mass in gamma (000) valley
constexpr double vel_const = (Efield * dirac)/(meff * m_e) ;

which is much clearer. The astute observer will notice that this fixes the first of my bugs: the mass of an electron is $9.11\times10^{-31}$ kg, not $9.11\times10^{-30}$ kg…

Dimensional Analysis

One way to verify scientific code is to check that the units of values used in computations are compatible. Making this possible with compile-time checks is the main motivation behind the Boost.Units library. For example, if I wanted to calculate the force on an electron in a given field:

constexpr auto echarge = si::constants::codata::e;
constexpr auto      cm = 0.01 * si::meter;

auto Efield = 1.2 * si::volt/cm;
quantity<si::force> f = Efield * echarge;

If I had the wrong units for electric field or electron charge, or I tried to store the result in a variable that was not a force quantity, I would get a compile error… as I did, when I refactored my separate energy calculation (second bug).

Keeping a Running Average

After each free flight period I need to update my running average velocity. I decided to replace my old manual averaging code with a statistical accumulator from the Boost.Accumulators library. It does mundane tasks like mean, but also allows you to extensively customize what you want calculated from your samples. At compile time it figures out the most efficient way to accomplish the various calculations you listed, simultaneously. For my simple case the usage is:

using namespace boost::accumulators;
accumulator_set<quantity<si::velocity>,  // samples are velocity type
                stats<tag::mean>>        // calculate only the mean
    vel_acc;

and on every sample:

vel_acc(cur_avg);                        // add avg velocity for this flight

Easy, right? Well, there is a little problem with it… Without modifications this produces some classic C++ “template spew”, the critical part of which is error: no match for ‘operator/’… It turns out that in order to calculate the mean, the Accumulators library divides the sum (of type quantity<si::velocity, double>) by the count (of type size_t), and the strongly typed Boost.Units does not permit this. As is often the case with Boost, the solution is to use a hook for adapting classes, in this case a specialization of boost::numeric::functional::fdiv, so the Units library understands the meaning of dividing these two types. Further details can be found in the Accumulators documentation.

Improving the Running Average

After some thought I realized that what I really needed was a weighted average velocity (bug 3). That is, I should consider the amount of time the electron was in free flight when I compute the mean. Fortunately with the Accumulators library that’s only a small change:

accumulator_set<quantity<si::velocity>,  // what we are storing
                stats<tag::mean>,        // what we want to calculate
                quantity<si::time>>      // the weight for each data point
    vel_acc;

and then for each sample we include the flight time:

vel_acc(cur_avg, weight = ts);

You may be unsurprised to hear that this change also required a specialization in the boost::numeric namespace - this time to use a floating-point $1.0$ as the default weight.

Conclusions

Though this is still fundamentally an academic exercise (many details and mechanisms of electron scattering were omitted or simplified), the process of modernizing the code was very satisfying. If you are interested in a more serious simulation, I found this Github repo which looks like a vastly more advanced version.

My code, warts and all, may be found here.

It was very enjoyable to modernize ancient code, and I recommend it. Maybe there is some in one of your projects!