Points: Measuring a Water Molecule

This tutorial introduces the experimental point component — Point, PointView, PointSet, and PointSetView, all in the chemist::experimental namespace — by working through a single problem:

Given the three nuclear positions of a water molecule, move the molecule so that its geometric center sits at the origin, measure its bond lengths and bond angle, and then hand its coordinates to an external library which knows nothing about Chemist.

Every code block below is compiled and run as part of Chemist’s test suite; see tests/cxx/doc_snippets/experimental_point_example.cpp and tests/python/doc_snippets/test_experimental_point_example.py.

Note

These classes live in chemist::experimental and are being developed alongside the existing chemist::Point<T>. Nothing in chemist:: has changed. The reasoning behind the design is written up in docs/source/developer/design/point/point.rst.

Getting a Number Back Out

Coordinates are stored type-erased, which means a point never has to expose whether it uses float or double, and code which takes a point never has to be templated on that choice. Similarly, asking a point for its x- coordinate gives you back a type-erased value. If you want an actual double you have to ask for it.

// The point classes deliberately keep the concrete floating-point type out of
// their API; a coordinate comes back as a type-erased value. Whichever type
// you actually want is yours to name, which is what this helper does.
double as_double(wtf::fp::FloatView<const wtf::fp::Float> value) {
    return value.value<double>();
}

Users of points are encouraged to avoid unwrapping the type-erased value unless they have to. Chemist is designed to work with Point objects directly to facilitate this.

One API, Two Ownership Models

A point can either own its coordinates or alias coordinates owned by something else. Point does the former and PointView does the latter, but they share a single API, so a function which computes a distance does not have to care which one it was handed:

C++Python

// Because Point and PointView share one API, a function which takes a
// read-only view works on both. A Point converts to one implicitly, so
// callers never write the conversion out.
double distance(const_point_view a, const_point_view b) {
    return as_double((a - b).magnitude());
}

const_point_view is an alias for PointView<const Point>, a view of a point you may read from, but not write to. Python has ImmutablePointView to serve the same purpose as const_point_view. Point converts to const_point_view (or in Python ImmutablePointView) implicitly, so callers pass a Point into functions expecting a const_point_view and the conversion happens on its own.

The same relationship holds one level up, between PointSet and PointSetView:

C++Python

// The same trick one level up: this takes a read-only view of a set, so it
// can be called with a PointSet, with a mutable view of one, or with a view
// over coordinates chemist does not own.
Point centroid_of(const_point_set_view points) {
    double x = 0.0, y = 0.0, z = 0.0;
    for(auto p : points) {
        x += as_double(p.get_x());
        y += as_double(p.get_y());
        z += as_double(p.get_z());
    }
    const auto n = static_cast<double>(points.size());
    return Point(x / n, y / n, z / n);
}

Building the Molecule

With those two helpers in hand, the molecule itself is unremarkable:

C++Python

    // A water molecule, in atomic units --- the same geometry the AO basis
    // set tutorial uses. The oxygen is first, then the two hydrogens.
    PointSet water{Point(0.0, -0.1432223429807816, 0.0),
                   Point(1.6380335020342418, 1.1365568803584036, 0.0),
                   Point(-1.6380335020342418, 1.1365568803584036, 0.0)};

    REQUIRE(water.size() == 3);

Internally, however, water is not a container of three Point objects. It holds three arrays — one for every point’s x-coordinate, one for every y, one for every z. That layout is what many HPC libraries expect because it vectorizes, and we will take advantage of it at the end of this tutorial.

That leaves the set with no Point objects inside it to hand out. Indexing into a PointSet gives you a PointView instead:

C++Python

    // The set stores each Cartesian direction in its own contiguous array,
    // but it behaves like a container of points. Indexing it gives you
    // something which acts exactly like a Point.
    auto oxygen = water[0];
    REQUIRE(as_double(oxygen.get_y()) == -0.1432223429807816);

    // What you got back is not a copy. It aliases the set, so writing through
    // it writes into the set.
    oxygen.set_z(1.0);
    REQUIRE(as_double(water[0].get_z()) == 1.0);
    oxygen.set_z(0.0);

The important part is the second half. Indexing the set does not copy a point out of it. You get a handle onto the set’s own storage, and writing through that handle writes into the set.

Centering the Molecule

That property is what makes translating the molecule a three-line loop. Each p is a view of the set, so each set_* call lands in the set’s coordinate arrays directly:

C++Python

    // Move the molecule so that its geometric center sits at the origin.
    // Every write goes straight into the set, because every `p` is a view of
    // it rather than a copy of one of its points.
    auto center = centroid_of(water);

    for(auto p : water) {
        p.set_x(as_double(p.get_x()) - as_double(center.get_x()));
        p.set_y(as_double(p.get_y()) - as_double(center.get_y()));
        p.set_z(as_double(p.get_z()) - as_double(center.get_z()));
    }

    // The centroid of the moved molecule is the origin.
    const auto off_center = as_double(centroid_of(water).magnitude());
    REQUIRE(off_center == Catch::Approx(0.0).margin(1.0e-15));

Note that centroid_of was written to take a read-only view of a set, and is being called here with a PointSet. As with distance, the conversion is implicit.

Measuring the Molecule

Bond lengths and bond angles come out of the same shared API. Subtracting two points gives the vector between them, and that vector is an owning Point even when both operands were views — so it stays valid after the expression which produced it:

C++Python

    // The two O-H bond lengths. distance() was written against a view, and
    // water[i] is one, so it can be called directly.
    const auto r1 = distance(water[1], water[0]);
    const auto r2 = distance(water[2], water[0]);
    REQUIRE(r1 == Catch::Approx(r2));
    REQUIRE(r1 == Catch::Approx(2.0786987791109155));

    // The H-O-H angle, from the inner product of the two bond vectors.
    // Subtracting two views gives an owning Point, so the vectors outlive the
    // expression which made them.
    const auto v1 = water[1] - water[0];
    const auto v2 = water[2] - water[0];

    const auto cos_theta =
      as_double(v1.inner_product(v2)) /
      (as_double(v1.magnitude()) * as_double(v2.magnitude()));
    // Water's H-O-H angle, a little under 104 degrees for this geometry.
    const auto degrees = std::acos(cos_theta) * 180.0 / std::acos(-1.0);
    REQUIRE(degrees == Catch::Approx(103.99968755694901));

Handing the Coordinates Out

Finally, the reason the set is laid out the way it is. Suppose an external library wants the coordinates as three bare arrays:

C++Python

// Stands in for an external library: it takes bare arrays of coordinates, one
// per Cartesian direction, and knows nothing about chemist.
double nuclear_repulsion(const double* xs, const double* ys, const double* zs,
                         const double* qs, std::size_t n) {
    double rv = 0.0;
    for(std::size_t i = 0; i < n; ++i) {
        for(std::size_t j = i + 1; j < n; ++j) {
            const auto dx = xs[i] - xs[j];
            const auto dy = ys[i] - ys[j];
            const auto dz = zs[i] - zs[j];
            rv += qs[i] * qs[j] / std::sqrt(dx * dx + dy * dy + dz * dz);
        }
    }
    return rv;
}

Because the set already stores one array per Cartesian direction, satisfying that request costs nothing. There is no repacking step; you ask the set for the array you want and take a pointer to it:

C++Python

    // An external library wants bare arrays. Because the set is stored as one
    // array per Cartesian direction, there is nothing to repack: ask it for
    // the array you want and hand over the pointer.
    auto xs = water.get_x_buffer();
    auto ys = water.get_y_buffer();
    auto zs = water.get_z_buffer();

    // The buffers are type-erased. Naming `double` here is the consumer's
    // obligation, not chemist's, and the buffer will say no if it is holding
    // something else. Contiguity is worth checking before taking a pointer.
    REQUIRE(xs.is_contiguous());

    auto x_span = xs.value<double>();
    auto y_span = ys.value<double>();
    auto z_span = zs.value<double>();

    // Nuclear charges, in the same order as the points.
    const double charges[] = {8.0, 1.0, 1.0};
    const auto v_nn        = nuclear_repulsion(x_span.data(), y_span.data(),
                                               z_span.data(), charges, water.size());

    // The nuclear repulsion energy of this geometry, in Hartree. Note that it
    // did not change when the molecule was translated.
    REQUIRE(v_nn == Catch::Approx(8.0023669741662));

Warning

A pointer obtained this way is only valid while the set’s storage is unchanged. Anything which grows the set — push_back, for instance — may reallocate and invalidate both the pointer and every PointView you are holding.

Passing Sets Around

The set-level views close the loop. A PointSetView is a handle onto somebody else’s points: it can read and write them, but it cannot grow or shrink the set, which makes it the right thing to hand to code that has no business resizing your molecule. The read-only flavor cannot even write:

C++Python

    // A view of the set is a handle, not a copy. Handing one out lets a caller
    // read or write the original without being able to grow or shrink it.
    point_set_view handle(water);
    handle[0].set_z(0.5);
    REQUIRE(as_double(water[0].get_z()) == 0.5);

    // A read-only view can not be written through at all; `set_z` does not
    // exist on the points it hands out, so the line below would not compile:
    //
    //     const_point_set_view(water)[0].set_z(0.5);
    const_point_set_view read_only(water);
    REQUIRE(read_only.size() == water.size());

    // When you really do want a copy, ask for one.
    auto copy = read_only.as_point_set();
    copy[0].set_z(99.0);
    REQUIRE(as_double(water[0].get_z()) == 0.5);

Recap

• Point owns three coordinates; PointView aliases three coordinates owned by something else. They share one API.

• PointSet owns three arrays of coordinates, one per Cartesian direction; PointSetView aliases them. They share one API.

• Indexing a set gives a view into it, not a copy out of it.

• Point converts implicitly to PointView, PointSet converts implicitly to PointSetView, and a mutable view converts implicitly to a read-only one. The reverse conversions do not exist.

• The concrete floating-point type is never part of the API. Consumers which need one name it themselves.