C++ Tutorials

This section present some hand-on tutorials to the chemfiles C++ API. All the code here is under the CC-0 Universal Licence which means that you are free to do whatever you want with it (i.e. it is Public Domain code)

Read a single frame

In this tutorials we will read a frame from a trajectory, and print the indexes of all the atom in the half-space x < 5.

We start by including the headers we need: chemfiles header, iostream for printing and vector for storing the indexes. All the public classes of chemfiles are accessible through the chemfiles.hpp header, which is the only header needed.

#include <chemfiles.hpp>
#include <iostream>
#include <vector>

Then we open a Trajectory and read the first frame:

chemfiles::Trajectory file("filename.xyz");
chemfiles::Frame frame = file.read();

We can now create a vector to store the indices of the atoms with x < 5, and get the positions of the atoms in the frame with the chemfiles::Frame::positions() function

std::vector<size_t> less_than_five;
auto positions = frame.positions();

Iterating through the atoms in the frame, we get the ones matching our condition. Frame::size() gives the number of atoms in the frame, which is also the size of the positions array. This array contains chemfiles::Vector3D representing the positions of the atoms in Angstroms.

for (size_t i=0; i<frame.size(); i++) {
    if (positions[i][0] < 5) {
        less_than_five.push_back(i);
    }
}

And finally we can print our results

std::cout << "Atoms with x < 5: " << std::endl;
for (auto i : less_than_five) {
    std::cout << "  - " << i << std::endl;
}

Click here to see the whole program

#include <chemfiles.hpp>
#include <iostream>
#include <vector>

int main() {
    chemfiles::Trajectory file("filename.xyz");
    chemfiles::Frame frame = file.read();

    std::vector<size_t> less_than_five;
    auto positions = frame.positions();

    for (size_t i=0; i<frame.size(); i++) {
        if (positions[i][0] < 5) {
            less_than_five.push_back(i);
        }
    }

    std::cout << "Atoms with x < 5: " << std::endl;
    for (auto i : less_than_five) {
        std::cout << "  - " << i << std::endl;
    }

    return 0;
}

For more information about reading frame in a trajectory, see the following functions:

• Trajectory::nsteps() and Trajectory::done() to know when to stop reading
• Trajectory::read_step() to directlty read a given step.
• Trajectory::set_cell() and Trajectory::set_topology() to specify an unit cell or a topology for all frames in a trajectory.

Generating a structure

Now that we know how to read frames from a trajectory, let’s try to create a new structure and write it to a file. As previsouly, we start by including the chemfiles header:

#include <chemfiles.hpp>

We start with a chemfiles::Frame that we construct by giving it the chemfiles::UnitCell that defines the periodic boundary conditions.

auto frame = chemfiles::Frame(chemfiles::UnitCell(10, 10, 10));

We can then add three chemfiles::Atom to this frame with their positions; and set the bonds between them to define a water molecule.

frame.add_atom(chemfiles::Atom("H"), {1, 0, 0});
frame.add_atom(chemfiles::Atom("O"), {0, 0, 0});
frame.add_atom(chemfiles::Atom("H"), {0, 1, 0});

frame.add_bond(0, 1);
frame.add_bond(2, 1);

Now that our frame is constructed, it is time to write it to a file. For that, we open a trajectory in write ('w') mode, and write to it.

auto trajectory = chemfiles::Trajectory("water.pdb", 'w');
trajectory.write(frame);

Click here to see the whole program

#include <chemfiles.hpp>

int main() {
    auto frame = chemfiles::Frame(chemfiles::UnitCell(10, 10, 10));

    frame.add_atom(chemfiles::Atom("H"), {1, 0, 0});
    frame.add_atom(chemfiles::Atom("O"), {0, 0, 0});
    frame.add_atom(chemfiles::Atom("H"), {0, 1, 0});

    frame.add_bond(0, 1);
    frame.add_bond(2, 1);

    auto trajectory = chemfiles::Trajectory("water.pdb", 'w');
    trajectory.write(frame);

    return 0;
}

Using selections

Now that we know how to read and write frame from trajectories, how about we do a bit a filtering? In this tutorial, we will read all the frames from a file, and use selections to filter which atoms we will write back to another file. This example will also show how chemfiles can be used to convert from a file format to another one.

We start by opening the two trajectories we will need

auto input = chemfiles::Trajectory("input.arc");
auto output = chemfiles::Trajectory("output.pdb", 'w');

And we create a chemfiles::Selection object to filter the atoms we want to keep.

auto selection = chemfiles::Selection("name Zn or name N");

Then we can iterate over all the frames in the trajectory

for (size_t step=0; step<input.nsteps(); step++) {
    auto frame = input.read();

And use the selection to get the list of atoms to remove. The result of Selection::list() is a std::vector<size_t> containing the atoms matching the selection.

auto to_remove = selection.list(frame);

In order to remove the atoms from the frame, we need to sort the to_remove vector in descending order: removing the atom at index i will shift the index of all the atoms after i. So we start from the end and work toward the start of the frame.

std::sort(std::begin(to_remove), std::end(to_remove), std::greater<size_t>());
for (auto i: to_remove) {
    frame.remove(i);
}

Finally, we can write the cleaned frame to the output file, and start the next iteration.

output.write(frame);

Click here to see the whole program

#include <chemfiles.hpp>
#include <functional>

int main() {
    auto input = chemfiles::Trajectory("input.arc");
    auto output = chemfiles::Trajectory("output.pdb", 'w');

    auto selection = chemfiles::Selection("name Zn or name N");

    for (size_t step=0; step<input.nsteps(); step++) {
        auto frame = input.read();
        auto to_remove = selection.list(frame);
        std::sort(std::begin(to_remove), std::end(to_remove), std::greater<size_t>());
        for (auto i: to_remove) {
            frame.remove(i);
        }
        output.write(frame);
    }

    return 0;
}