Atomic units

This example is different from the others. Instead of showing how Au can improve raw numeric code, this one shows how existing users of Au can extend it to support a new system of units.

Most units — including all units shipped with Au — are defined by exact relationships with the SI base units. The Hartree atomic units in this example are very different, in an interesting way:

• Within Hartree atomic units, relationships are (still) exact.

• Between Hartree atomic units and SI units, relationships are relative to measured constants.

These constants have associated measurement uncertainty, and future improvements in physics will surely update their values. Unfortunately, Au’s entire machinery is built on representing exact relationships between units, not approximate or uncertain ones.

This makes it sound like Hartree units should be impossible with Au. However, we’ll see that in fact, we can handle them just fine! The one thing we can’t do is ship these definitions with the library. That would mean choosing one-size-fits-all values for these uncertain and measured constants, preventing each project from choosing the measurement version best suited to its needs.

Tip

This example provides a complete demonstration of creating a new system of units. You can use it for inspiration if your project has any units with inexact relationships. And if your project uses Hartree atomic units specifically, you can even copy paste these files and use them directly!

The definitions

Hartree atomic units are the working units of quantum chemistry. They are convenient precisely because they are defined so that the quantities in the Schrödinger equation come out around 1.

The definitions themselves are a library of two files — a header (.hh) and an implementation file (.cc) — which a program then includes. We show all three below.

Note

If you are using C++17 or later, you can skip the .cc file entirely, by using inline variable definitions for the labels. See How to define new units for details.

Header (atomic_units.hh)

We’ll first present the code, and then discuss the design choices that went into it.

namespace atomic_units {

//
// Exact inputs: SI-defining constants.  These have no uncertainty; the SI *defines* them.
//

constexpr auto c = au::SPEED_OF_LIGHT;
constexpr auto hbar = au::REDUCED_PLANCK_CONSTANT;
constexpr auto e = au::ELEMENTARY_CHARGE;

namespace detail {

//
// Measured inputs.  These are the only two values here that a future experiment can revise.
//

// CODATA 2018: 7.297 352 5693(11) e-3
constexpr auto fine_structure_constant() {
    using namespace ::au::au_literals;
    return 7.297'352'5693e-3_mag;
}

// CODATA 2018: 9.109 383 7015(28) e-31 kg
constexpr auto electron_mass() {
    using namespace ::au::au_literals;
    using au::kilo;
    return kilo(9.109'383'7015e-31_g);
}

}  // namespace detail

// Fine structure constant.
constexpr auto alpha = detail::fine_structure_constant();

// Electron mass.  Relabeled `m_e`: `associated_unit` recovers the unit the literal named, and
// `make_constant` below rebuilds the constant on top of it.
struct ElectronMasses : decltype(associated_unit(detail::electron_mass())) {
    static constexpr const char label[] = "m_e";
};
constexpr auto ELECTRON_MASS = au::make_constant(ElectronMasses{});
constexpr auto m_e = ELECTRON_MASS;

//
// Derived atomic units.  Each is defined by its physical formula, not by a decimal value, so the
// relationships among them are exact.
//

// Energy: the hartree, E_h = m_e c^2 alpha^2.
struct Hartrees : decltype(associated_unit(m_e * squared(c * alpha))) {
    static constexpr const char label[] = "E_h";
};
constexpr auto HARTREE = au::make_constant(Hartrees{});
constexpr auto hartree = au::SingularNameFor<Hartrees>{};
constexpr auto hartrees = au::QuantityMaker<Hartrees>{};

// Length: the Bohr radius, a_0 = hbar / (m_e c alpha).
struct BohrRadii : decltype(associated_unit(hbar / (m_e * c * alpha))) {
    static constexpr const char label[] = "a_0";
};
constexpr auto BOHR_RADIUS = au::make_constant(BohrRadii{});
constexpr auto bohr_radius = au::SingularNameFor<BohrRadii>{};
constexpr auto bohr_radii = au::QuantityMaker<BohrRadii>{};

// Time: t_a = hbar / E_h.
struct AtomicTimeUnits : decltype(associated_unit(hbar / hartrees)) {
    static constexpr const char label[] = "t_a";
};
constexpr auto ATOMIC_TIME = au::make_constant(AtomicTimeUnits{});
constexpr auto atomic_time_unit = au::SingularNameFor<AtomicTimeUnits>{};
constexpr auto atomic_time_units = au::QuantityMaker<AtomicTimeUnits>{};

// Charge: the elementary charge itself.  This one needs no `struct` of its own: `e` is already
// exact, and Au's unit for it already carries the label `e`, so there is nothing to derive and
// nothing to relabel.  We only name the makers, to round out the system.
using ElementaryCharges = decltype(associated_unit(e));
constexpr auto elementary_charge = au::SingularNameFor<ElementaryCharges>{};
constexpr auto elementary_charges = au::QuantityMaker<ElementaryCharges>{};

//
// This `static_assert` proves that our unit definitions are exact, within the system.  If it were
// otherwise, this conversion to `int` would not compile.
//
static_assert(hbar.in<int>(hartrees * atomic_time_units) == 1,
              "One reduced Planck constant must be exactly one hartree times one atomic time unit");

}  // namespace atomic_units

Since this is a header file, we qualify Au names with the au:: prefix, rather than importing them. This avoids namespace pollution; see Namespaces and includes for more details. (The one exception is utilities found by argument-dependent lookup, discussed below.)

The header also needs the two measured values. By far the best way to express them is using Au literals: the magnitude literal 7.297'352'5693e-3_mag, and the prefixed unit literal kilo(9.109'383'7015e-31_g). This lets the source code carry the published CODATA digits — the internationally recommended values of the fundamental physical constants, revised every few years — verbatim, making it easy to check at a glance.¹ Since it’s a header file, we do still need to avoid namespace pollution, so each one goes inside its own function, where the using-directive expires at the closing brace and never reaches a consumer.

From this point, we simply follow our standard approach for defining new units and constants. One key thing to note is that when we define a new unit, the type it inherits from must be a unit, not a constant or anything else. The most ergonomic way to do this is to write the most convenient expression, and then pass it through the associated_unit utility: this converts any unit slot into an actual unit type.

Also, as a minor point of style, note that we do not need the au:: prefix on utilities that take au:: types as arguments (associated_unit, squared). The compiler uses argument-dependent lookup (ADL) to find these functions.

Implementation (atomic_units.cc)

In C++14, the unit labels need one out-of-line definition each, in a .cc file:

// C++14 needs these out-of-line definitions for the unit labels.  In C++17 or later, declare the
// labels `static constexpr inline` in the header instead, and delete this file.
namespace atomic_units {

constexpr const char ElectronMasses::label[];
constexpr const char Hartrees::label[];
constexpr const char BohrRadii::label[];
constexpr const char AtomicTimeUnits::label[];

}  // namespace atomic_units

Using them

Here’s an example showing how you would use these definitions in your program. As usual, we provide “front matter” (includes and using statements) collapsed by default, and then show the core code itself.

Includes and usings

#include <iostream>

#include "au/io.hh"
#include "au/units/coulombs.hh"
#include "au/units/joules.hh"
#include "au/units/meters.hh"
#include "au/units/seconds.hh"
#include "examples/atomic_units/atomic_units.hh"

// This is a `.cc` file, so we import the names we use, one at a time.  Note that the header this
// includes does the opposite: it qualifies Au names with `au::` rather than importing any, because
// a namespace-scope `using` in a header would leak into every file that includes it.  See the
// "Namespaces and includes" discussion page.
using au::coulombs;
using au::joules;
using au::meters;
using au::seconds;

using atomic_units::ATOMIC_TIME;
using atomic_units::atomic_time_units;
using atomic_units::BOHR_RADIUS;
using atomic_units::e;
using atomic_units::HARTREE;
using atomic_units::hartrees;
using atomic_units::hbar;

int main() {
    // Crossing out to SI.  The accuracy here is set by the two measured inputs, and nothing else:
    // Au composes the entire definition chain into one exact rational factor before applying it.
    //
    // Note that no unit label is written out below.  Every name in this output --- `a_0`, `m`, and
    // the rest --- is the label Au derives from the unit itself.
    std::cout << BOHR_RADIUS << " = " << BOHR_RADIUS.as<double>(meters) << '\n';
    std::cout << HARTREE << " = " << HARTREE.as<double>(joules) << '\n';
    std::cout << ATOMIC_TIME << " = " << ATOMIC_TIME.as<double>(seconds) << '\n';

    // The atomic unit of charge is the one that crosses to SI *exactly*: the elementary charge is
    // an SI-defining constant, so no measured quantity enters this conversion at all.
    std::cout << e << " = " << e.as<double>(coulombs) << '\n';

    // Staying inside the system.  This prints `1`, not `0.9999999997`, because the atomic time unit
    // is *defined* as hbar / E_h rather than pasted in as a decimal.  (This claim also has stronger
    // evidence than just the printed output: the `static_assert` in the header file.)
    std::cout << hbar << " = " << hbar.as<double>(hartrees * atomic_time_units) << '\n';
}

This prints:

a_0 = 5.29177e-11 m
E_h = 4.35974e-18 J
t_a = 2.41888e-17 s
e = 1.60218e-19 C
h_bar = 1 E_h * t_a

Every name in that output is derived, not typed. a_0, E_h, t_a, e and h_bar come from the constants themselves; m, J, s, C and E_h * t_a come from the units they were converted into.

Note that the charge line is the only one whose conversion to SI is exact: the elementary charge is an SI-defining constant, so neither measured input is involved.

The last line is the interesting one, and it makes a different kind of claim. h_bar comes out as exactly 1 in units of E_h * t_a — a relationship inside the system, with no measured input anywhere in it. Below, we’ll see how the program gets that checked by the compiler, rather than asking you to trust a printed digit.

What’s happening

Let’s highlight the most salient points about how this new system works.

Everything derives from exactly five inputs. Three of them — the speed of light, the reduced Planck constant, and the elementary charge — are SI-defining constants. They have no uncertainty, because the SI defines the kilogram, meter and second in terms of them rather than the other way round. The other two — the electron mass and the fine structure constant — are measured. Separating them this way means that updating to a newer CODATA release is a two-line change, and you can see at a glance exactly which numbers in your program are experimental inputs.

Each unit is defined by its formula, not by a decimal value. The hartree is written as m_e * squared(c * alpha), not as 4.3597447e-18 joules. That is the difference that makes the last line of the output read 1 E_h * t_a rather than 0.999999999.... The atomic time unit is defined as $\hbar/E_h$, so within the system that relationship is exact — Au carries it as an exact rational magnitude, never as a rounded double.

The two core properties we outlined at the beginning follow from these definitions:

• Inside the system, relationships are exact. No accumulated error from round-tripping through SI, no drift, no epsilon comparisons.

• Crossing out to SI is as accurate as physics currently allows. Au composes the entire definition chain — through the hartree, through $\alpha^2$, through the electron mass — into a single exact rational factor, and applies it once. The only error is the uncertainty in the two measured inputs.

The first property is checked by the compiler: The header closes with a static_assert that asks for $\hbar$’s value in hartrees * atomic_time_units as an int.

static_assert(hbar.in<int>(hartrees * atomic_time_units) == 1,
              "One reduced Planck constant must be exactly one hartree times one atomic time unit");

The int output is what gives this teeth. Au refuses to compile a conversion it knows would truncate, so the fact that this line compiles at all means the conversion factor is exactly 1, and not a double that happens to round to 1.0. And this property is completely independent of the actual measured values, so it won’t break when we update them after newer experiments give us more precise values.

The labels make output readable. Because each unit declares a label, printing a quantity gives a_0 or E_h rather than an unwieldy composite of SI base units — and composite units get composite labels, which is where the E_h * t_a on the last line comes from. Streaming a Constant prints its label alone, so std::cout << HARTREE writes E_h. That is why the program above never spells a unit name out in a string literal.

Nothing here is special-cased. These are ordinary Au units, using the same mechanism as any other custom unit. They compose with the built-in units, participate in the same conversion checks, and cost nothing at runtime.

Related reading

• How to define new units, including the C++14 versus C++17 label forms.
• How to define new constants.
• Magnitude, the exact compile-time arithmetic that makes the in-system relationships exact.
• Unit literals and prefixed unit literals, for writing a measured value exactly as it was published.
• Applying a prefix to a Constant, which is what makes kilo(...e-31_g) mean what you would want it to mean.
• Namespaces and includes, including why a .cc imports names individually and a .hh does not.

---

1. Compare this readability to the alternative, mag<72973525693>() * pow<-13>(mag<10>()), with its manually shifted powers of 10. ↩