C++ Tutorials¶

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;
}

#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>

Everything starts in a chemfiles::Topology. This is the class that defines the atoms and the connectivity in a system. Here, we add three chemfiles::Atom and two bonds to create a water molecule.

chemfiles::Topology topology;
topology.add_atom(chemfiles::Atom("H"));
topology.add_atom(chemfiles::Atom("O"));
topology.add_atom(chemfiles::Atom("H"));

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

We can then create a chemfiles::Frame corresponding to this chemfiles::Topology and set the atomic positions.

chemfiles::Frame frame(topology);
auto positions = frame.positions();

positions[0] = chemfiles::Vector3D(1, 0, 0);
positions[1] = chemfiles::Vector3D(0, 0, 0);
positions[2] = chemfiles::Vector3D(0, 1, 0);

Another possibility is to directly add atoms to the frame. Here we define a second molecule representing carbon dioxyde. chemfiles::Frame::add_atom() takes two arguments: the atom, and the position of the atom as a chemfiles::Vector3D.

frame.add_atom(chemfiles::Atom("O"), {5, 0, 0});
frame.add_atom(chemfiles::Atom("C"), {6, 0, 0});
frame.add_atom(chemfiles::Atom("O"), {7, 0, 0});
frame.add_bond(3, 4);
frame.add_bond(4, 5);

Finally, we can set the chemfiles::UnitCell associated with this frame.

frame.set_cell(chemfiles::UnitCell(10, 10, 10));

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-co2.pdb", 'w');
trajectory.write(frame);

#include <chemfiles.hpp>

int main() {
    chemfiles::Topology topology;
    topology.add_atom(chemfiles::Atom("H"));
    topology.add_atom(chemfiles::Atom("O"));
    topology.add_atom(chemfiles::Atom("H"));

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

    chemfiles::Frame frame(topology);
    auto positions = frame.positions();

    positions[0] = chemfiles::Vector3D(1, 0, 0);
    positions[1] = chemfiles::Vector3D(0, 0, 0);
    positions[2] = chemfiles::Vector3D(0, 1, 0);

    frame.add_atom(chemfiles::Atom("O"), {5, 0, 0});
    frame.add_atom(chemfiles::Atom("C"), {6, 0, 0});
    frame.add_atom(chemfiles::Atom("O"), {7, 0, 0});
    frame.add_bond(3, 4);
    frame.add_bond(4, 5);

    frame.set_cell(chemfiles::UnitCell(10, 10, 10));

    auto trajectory = chemfiles::Trajectory("water-co2.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);

#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;
}