Frame

class chemfiles::Frame

A frame contains data from one simulation step The Frame class holds data from one step of a simulation: the current topology, the positions, and the velocities of the particles in the system. If some information is missing the corresponding data is filled with a default value. Specifically:

  • cell is an infinite unit cell;

  • topology is empty, and contains no data;

  • positions is filled with zeros;

  • velocities is the nullopt variant of chemfiles::optional.

Iterating over a Frame will yield all the atoms in the system.

auto frame = Frame();
frame.add_atom(Atom("Fe"), {0.0, 0.0, 0.0});
frame.add_atom(Atom("Fe"), {1.0, 1.0, 1.0});
frame.add_atom(Atom("Fe"), {2.0, 2.0, 2.0});

for (Atom& atom: frame) {
    assert(atom.name() == "Fe");
}

Public Functions

Frame(UnitCell cell = UnitCell())

Create an empty frame with no atoms and the given cell.

auto frame = Frame();
assert(frame.size() == 0);
assert(frame.cell() == UnitCell());

Frame clone() const

Get a clone (exact copy) of this frame.

This replace the implicit copy constructor (which is private) to make an explicit copy of the frame.

auto frame = Frame();
assert(frame.size() == 0);

auto copy = frame.clone();
copy.resize(42);

assert(frame.size() == 0);
assert(copy.size() == 42);

const Topology &topology() const

Get a const reference to the topology of this frame

It is not possible to get a modifiable reference to the topology, because it would then be possible to remove/add atoms without changing the actual positions and velocity storage. Instead, all the mutating functionalities of the topology are mirrored on the frame (adding and removing bonds, adding residues, etc.)

auto frame = Frame();
frame.resize(3);

auto topology = frame.topology();
assert(topology[0].name() == "");

// Manually constructing a topology
topology = Topology();
topology.add_atom(Atom("H"));
topology.add_atom(Atom("O"));
topology.add_atom(Atom("H"));
topology.add_bond(0, 1);
topology.add_bond(2, 1);

frame.set_topology(topology);
assert(frame.topology()[0].name() == "H");

void set_topology(Topology topology)

Set the topology of this frame to topology

auto frame = Frame();
frame.resize(3);

auto topology = frame.topology();
assert(topology[0].name() == "");

// Manually constructing a topology
topology = Topology();
topology.add_atom(Atom("H"));
topology.add_atom(Atom("O"));
topology.add_atom(Atom("H"));
topology.add_bond(0, 1);
topology.add_bond(2, 1);

frame.set_topology(topology);
assert(frame.topology()[0].name() == "H");

Parameters

  • topology: the new Topology to use for this frame

Exceptions

  • Error: if the topology size does not match the size of this frame

const UnitCell &cell() const

Get a const reference to the unit cell of this frame

auto frame = Frame();

auto cell = frame.cell();
assert(cell.shape() == UnitCell::INFINITE);

cell = UnitCell({23, 34, 11.5});
frame.set_cell(cell);

assert(frame.cell().shape() == UnitCell::ORTHORHOMBIC);

UnitCell &cell()

Get a reference to the unit cell of this frame

auto frame = Frame();

auto cell = frame.cell();
assert(cell.shape() == UnitCell::INFINITE);

cell = UnitCell({23, 34, 11.5});
frame.set_cell(cell);

assert(frame.cell().shape() == UnitCell::ORTHORHOMBIC);

void set_cell(UnitCell cell)

Set the unit cell for this frame to cell

auto frame = Frame();

auto cell = frame.cell();
assert(cell.shape() == UnitCell::INFINITE);

cell = UnitCell({23, 34, 11.5});
frame.set_cell(cell);

assert(frame.cell().shape() == UnitCell::ORTHORHOMBIC);

Parameters

  • cell: the new UnitCell to use for this frame

size_t size() const

Get the number of atoms in this frame

auto frame = Frame();
assert(frame.size() == 0);

frame.resize(10);
assert(frame.size() == 10);

span<Vector3D> positions()

Get the positions of the atoms in this frame.

auto frame = Frame();
frame.add_atom(Atom("H"), {3.0, 4.0, 5.0});
frame.add_atom(Atom("O"), {1.0, -2.0, 3.0});
frame.add_atom(Atom("H"), {1.3, 0.0, -1.0});

auto positions = frame.positions();
assert(positions.size() == 3);

// Indexing the positions
assert(positions[0] == Vector3D(3.0, 4.0, 5.0));
assert(positions[1] == Vector3D(1.0, -2.0, 3.0));
assert(positions[2] == Vector3D(1.3, 0.0, -1.0));

// Iteration on positions
for (auto& position: positions) {
    position[0] += 1.0;
    position[2] -= 1.0;
}

assert(positions[0] == Vector3D(4.0, 4.0, 4.0));
assert(positions[1] == Vector3D(2.0, -2.0, 2.0));
assert(positions[2] == Vector3D(2.3, 0.0, -2.0));

const std::vector<Vector3D> &positions() const

Get the positions in this frame as a const reference

auto frame = Frame();
frame.add_atom(Atom("H"), {3.0, 4.0, 5.0});
frame.add_atom(Atom("O"), {1.0, -2.0, 3.0});
frame.add_atom(Atom("H"), {1.3, 0.0, -1.0});

auto positions = frame.positions();
assert(positions.size() == 3);

// Indexing the positions
assert(positions[0] == Vector3D(3.0, 4.0, 5.0));
assert(positions[1] == Vector3D(1.0, -2.0, 3.0));
assert(positions[2] == Vector3D(1.3, 0.0, -1.0));

// Iteration on positions
for (auto& position: positions) {
    position[0] += 1.0;
    position[2] -= 1.0;
}

assert(positions[0] == Vector3D(4.0, 4.0, 4.0));
assert(positions[1] == Vector3D(2.0, -2.0, 2.0));
assert(positions[2] == Vector3D(2.3, 0.0, -2.0));

void add_velocities()

Add velocities data storage to this frame.

If velocities are already defined, this functions does nothing. The new velocities are initialized to 0.

auto frame = Frame();
// Default constructed frames does not contains velocities
assert(!frame.velocities());
frame.add_velocities();

assert(frame.velocities());

optional<span<Vector3D>> velocities()

Get an velocities of the atoms in this frame, if this frame contains velocity data.

auto frame = Frame();
// Default constructed frames does not contains velocities
assert(!frame.velocities());

frame.add_velocities();

// adding a few atoms with velocity data
frame.add_atom(Atom("H"), {}, {3.0, 4.0, 5.0});
frame.add_atom(Atom("O"), {}, {1.0, -2.0, 3.0});
frame.add_atom(Atom("H"), {}, {1.3, 0.0, -1.0});

auto velocities = *frame.velocities();
assert(velocities.size() == 3);

// Indexing the velocities
assert(velocities[0] == Vector3D(3.0, 4.0, 5.0));
assert(velocities[1] == Vector3D(1.0, -2.0, 3.0));
assert(velocities[2] == Vector3D(1.3, 0.0, -1.0));

// Iteration on positions
for (auto& velocity: velocities) {
    velocity[0] += 1.0;
    velocity[2] -= 1.0;
}

assert(velocities[0] == Vector3D(4.0, 4.0, 4.0));
assert(velocities[1] == Vector3D(2.0, -2.0, 2.0));
assert(velocities[2] == Vector3D(2.3, 0.0, -2.0));

optional<const std::vector<Vector3D>&> velocities() const

Get an velocities of the atoms in this frame as a const reference, if this frame contains velocity data.

auto frame = Frame();
// Default constructed frames does not contains velocities
assert(!frame.velocities());

frame.add_velocities();

// adding a few atoms with velocity data
frame.add_atom(Atom("H"), {}, {3.0, 4.0, 5.0});
frame.add_atom(Atom("O"), {}, {1.0, -2.0, 3.0});
frame.add_atom(Atom("H"), {}, {1.3, 0.0, -1.0});

auto velocities = *frame.velocities();
assert(velocities.size() == 3);

// Indexing the velocities
assert(velocities[0] == Vector3D(3.0, 4.0, 5.0));
assert(velocities[1] == Vector3D(1.0, -2.0, 3.0));
assert(velocities[2] == Vector3D(1.3, 0.0, -1.0));

// Iteration on positions
for (auto& velocity: velocities) {
    velocity[0] += 1.0;
    velocity[2] -= 1.0;
}

assert(velocities[0] == Vector3D(4.0, 4.0, 4.0));
assert(velocities[1] == Vector3D(2.0, -2.0, 2.0));
assert(velocities[2] == Vector3D(2.3, 0.0, -2.0));

void resize(size_t size)

Resize the frame to contain size atoms.

If the new number of atoms is bigger than the old one, missing data is initialized to 0. Pre-existing values are conserved.

If the new size if smaller than the old one, all atoms and connectivity elements after the new size are removed.

auto frame = Frame();
frame.resize(10);
assert(frame.size() == 10);

// new atoms contains default data
for (auto position: frame.positions()) {
    assert(position[0] == 0.0);
    assert(position[1] == 0.0);
    assert(position[2] == 0.0);
}

for (auto atom: frame.topology()) {
    assert(atom.name() == "");
}

void reserve(size_t size)

Allocate memory in the frame to have enough size for size atoms.

This function does not change the actual number of atoms in the frame, and should be used as an optimisation.

auto frame = Frame();
frame.resize(10);
assert(frame.size() == 10);

// reserve allocate memory, but does not change the size
frame.reserve(100);
assert(frame.size() == 10);

void add_atom(Atom atom, Vector3D position, Vector3D velocity = Vector3D())

Add an atom at the given position and optionally with the given velocity. The velocity value will only be used if this frame contains velocity data.

auto frame = Frame();
// add atom without velocities
frame.add_atom(Atom("H"), Vector3D(3.0, 4.0, 5.0));

frame.add_velocities();
// add atom with velocities
frame.add_atom(Atom("O"), {0.0, 0.0, 0.0}, {1.0, 2.0, 0.0});

auto velocities = *frame.velocities();
assert(velocities[0] == Vector3D(0.0, 0.0, 0.0));
assert(velocities[1] == Vector3D(1.0, 2.0, 0.0));

// using brace initialization for smaller code
frame.add_atom(Atom("H"), {1.0, 2.0, 0.0}, {0.0, -0.4, 0.3});

void remove(size_t i)

Remove the atom at index i in the system.

auto frame = Frame();
frame.add_atom(Atom("H"), {1.0, 0.0, 0.0});
frame.add_atom(Atom("O"), {0.0, 1.0, 0.0});
frame.add_atom(Atom("H"), {0.0, 0.0, 1.0});
assert(frame.size() == 3);

assert(frame.topology()[1].name() == "O");
assert(frame.positions()[1] == Vector3D(0.0, 1.0, 0.0));

frame.remove(1);
assert(frame.size() == 2);

// Removing an atom changes the indexes of atoms after the one removed
assert(frame.topology()[1].name() == "H");
assert(frame.positions()[1] == Vector3D(0.0, 0.0, 1.0));

Exceptions

size_t step() const

Get the current simulation step.

The step is set by the Trajectory when reading a frame.

auto frame = Frame();
assert(frame.step() == 0);

frame.set_step(424);
assert(frame.step() == 424);

void set_step(size_t step)

Set the current simulation step to step

auto frame = Frame();
assert(frame.step() == 0);

frame.set_step(424);
assert(frame.step() == 424);

void guess_bonds()

Guess the bonds, angles, dihedrals and impropers angles in this frame.

The bonds are guessed using a distance-based algorithm, and then angles, dihedrals and impropers are guessed from the bonds. The distance criterion uses the Van der Waals radii of the atoms. If this information is missing for a specific atoms, one can use configuration files to provide it.

// Building a frame containing a Cl2 molecule
auto frame = Frame();
frame.add_atom(Atom("Cl"), {0.0, 0.0, 0.0});
frame.add_atom(Atom("Cl"), {2.0, 0.0, 0.0});

assert(frame.topology().bonds().size() == 0);

frame.guess_bonds();
assert(frame.topology().bonds().size() == 1);

Exceptions

  • Error: if the Van der Waals radius in unknown for a given atom.

void clear_bonds()

Remove all connectivity information in the frame’s topology

auto frame = Frame();
frame.add_atom(Atom("H"), {1.0, 0.0, 0.0});
frame.add_atom(Atom("O"), {0.0, 0.0, 0.0});
frame.add_atom(Atom("H"), {0.0, 1.0, 0.0});

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

assert(frame.topology().bonds().size() == 2);
assert(frame.topology().angles().size() == 1);

frame.clear_bonds();
assert(frame.topology().bonds().size() == 0);
assert(frame.topology().angles().size() == 0);

void add_residue(Residue residue)

Add a residue to this frame’s topology.

auto frame = Frame();
frame.add_atom(Atom("Zn"), {0.0, 0.0, 0.0});
frame.add_atom(Atom("Fe"), {1.0, 2.0, 3.0});

auto residue = Residue("first");
residue.add_atom(0);
frame.add_residue(residue);

// residues are actually stored in the topology
assert(frame.topology().residues().size() == 1);

Parameters

  • residue: the residue to add to this topology

Exceptions

  • chemfiles::Error: if any atom in the residue is already in another residue in this topology. In that case, the topology is not modified.

void add_bond(size_t atom_i, size_t atom_j, Bond::BondOrder bond_order = Bond::UNKNOWN)

Add a bond in the system, between the atoms at index atom_i and atom_j.

auto frame = Frame();
frame.add_atom(Atom("H"), {1.0, 0.0, 0.0});
frame.add_atom(Atom("O"), {0.0, 0.0, 0.0});
frame.add_atom(Atom("H"), {0.0, 1.0, 0.0});

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

// the bonds are actually stored inside the topology
assert(frame.topology().bonds() == std::vector<Bond>({{0, 1}, {1, 2}}));
// angles are automaticaly computed too
assert(frame.topology().angles() == std::vector<Angle>({{0, 1, 2}}));

Parameters

  • atom_i: the index of the first atom in the bond

  • atom_j: the index of the second atom in the bond

  • bond_order: the bond order of the new bond

Exceptions

  • OutOfBounds: if atom_i or atom_j are greater than size()

  • Error: if atom_i == atom_j, as this is an invalid bond

void remove_bond(size_t atom_i, size_t atom_j)

Remove a bond in the system, between the atoms at index atom_i and atom_j.

If the bond does not exist, this does nothing.

auto frame = Frame();
frame.add_atom(Atom("H"), {1.0, 0.0, 0.0});
frame.add_atom(Atom("O"), {0.0, 0.0, 0.0});
frame.add_atom(Atom("H"), {0.0, 1.0, 0.0});

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

assert(frame.topology().bonds() == std::vector<Bond>({{0, 1}, {1, 2}}));

frame.remove_bond(1, 0);
assert(frame.topology().bonds() == std::vector<Bond>({{1, 2}}));

// This does nothing
frame.remove_bond(0, 2);
assert(frame.topology().bonds() == std::vector<Bond>({{1, 2}}));

Parameters

  • atom_i: the index of the first atom in the bond

  • atom_j: the index of the second atom in the bond

Exceptions

Atom &operator[](size_t index)

Get a reference to the atom at the position index.

auto frame = Frame();
frame.add_atom(Atom("Co"), {0.0, 0.0, 0.0});
frame.add_atom(Atom("V"), {1.0, 0.0, 0.0});
frame.add_atom(Atom("Fe"), {0.0, 2.0, 0.0});
frame.add_atom(Atom("Fe"), {0.0, 0.0, 3.0});

assert(frame[0].name() == "Co");
assert(frame[1].name() == "V");

frame[2].set_mass(45);
assert(frame[2].mass() == 45);

Parameters

  • index: the atomic index

Exceptions

const Atom &operator[](size_t index) const

Get a const reference to the atom at the position index.

auto frame = Frame();
frame.add_atom(Atom("Co"), {0.0, 0.0, 0.0});
frame.add_atom(Atom("V"), {1.0, 0.0, 0.0});
frame.add_atom(Atom("Fe"), {0.0, 2.0, 0.0});
frame.add_atom(Atom("Fe"), {0.0, 0.0, 3.0});

assert(frame[0].name() == "Co");
assert(frame[1].name() == "V");

frame[2].set_mass(45);
assert(frame[2].mass() == 45);

Parameters

  • index: the atomic index

Exceptions

double distance(size_t i, size_t j) const

Get the distance between the atoms at indexes i and j, accounting for periodic boundary conditions. The distance is expressed in angstroms.

auto frame = Frame();
frame.add_atom(Atom(""), {0.0, 0.0, 0.0});
frame.add_atom(Atom(""), {1.0, 2.0, 3.0});

assert(fabs(frame.distance(0, 1) - sqrt(14)) < 1e-15);

Exceptions

double angle(size_t i, size_t j, size_t k) const

Get the angle formed by the atoms at indexes i, j and k, accounting for periodic boundary conditions. The angle is expressed in radians.

auto frame = Frame();
frame.add_atom(Atom(""), {1.0, 0.0, 0.0});
frame.add_atom(Atom(""), {0.0, 0.0, 0.0});
frame.add_atom(Atom(""), {0.0, 1.0, 0.0});

assert(fabs(frame.angle(0, 1, 2) - M_PI / 2) < 1e-12);

Exceptions

double dihedral(size_t i, size_t j, size_t k, size_t m) const

Get the dihedral angle formed by the atoms at indexes i, j, k and m, accounting for periodic boundary conditions. The angle is expressed in radians.

auto frame = Frame();
frame.add_atom(Atom(""), {1.0, 0.0, 0.0});
frame.add_atom(Atom(""), {0.0, 0.0, 0.0});
frame.add_atom(Atom(""), {0.0, 1.0, 0.0});
frame.add_atom(Atom(""), {0.0, 1.0, 1.0});

assert(fabs(frame.dihedral(0, 1, 2, 3) - M_PI / 2) < 1e-12);

Exceptions

double out_of_plane(size_t i, size_t j, size_t k, size_t m) const

Get the out of plane distance formed by the atoms at indexes i, j, k and m, accounting for periodic boundary conditions. The distance is expressed in angstroms.

This is the distance between the atom j and the ikm plane. The j atom is the center of the improper dihedral angle formed by i, j, k and m.

auto frame = Frame();
frame.add_atom(Atom(""), {0.0, 0.0, 0.0});
frame.add_atom(Atom(""), {0.0, 0.0, 2.0});
frame.add_atom(Atom(""), {1.0, 0.0, 0.0});
frame.add_atom(Atom(""), {0.0, 1.0, 0.0});

assert(frame.out_of_plane(0, 1, 2, 3) == 2);

Exceptions

const property_map &properties() const

Get the map of properties associated with this frame. This map might be iterated over to list the properties of the frame, or directly accessed.

auto frame = Frame();

frame.set("a string", "the lazy fox");
frame.set("a number", 122);

// Iterator over properties in the frame
for (auto it: frame.properties()) {
    if (it.first == "a string") {
        assert(it.second.as_string() == "the lazy fox");
    } else if (it.first == "a number") {
        assert(it.second.as_double() == 122);
    }
}

void set(std::string name, Property value)

Set an arbitrary property for this frame with the given name and value. If a property with this name already exist, it is silently replaced with the new value.

auto frame = Frame();

frame.set("foo", Property(-23));
assert(frame.get("foo")->as_double() == -23);

// Override the 'foo' property
frame.set("foo", Property(false));
assert(frame.get("foo")->as_bool() == false);

optional<const Property&> get(const std::string &name) const

Get the Property with the given name for this frame if it exists.

If no property with the given name is found, this function returns nullopt.

auto frame = Frame();
frame.set("foo", Property(23));

assert(frame.get("foo"));
assert(frame.get("foo")->as_double() == 23);

assert(frame.get<Property::DOUBLE>("foo").value() == 23);
assert(!frame.get<Property::STRING>("foo"));

assert(!frame.get("bar"));

template<Property::Kind kind>
optional<typename property_metadata<kind>::type> get(const std::string &name) const

Get the Property with the given name for this frame if it exists, and check that it has the required kind.

If no property with the given name is found, this function returns nullopt.

If a property with the given name is found, but has a different kind, this function emits a warning and returns nullopt.

auto frame = Frame();
frame.set("foo", Property(23));

assert(frame.get("foo"));
assert(frame.get("foo")->as_double() == 23);

assert(frame.get<Property::DOUBLE>("foo").value() == 23);
assert(!frame.get<Property::STRING>("foo"));

assert(!frame.get("bar"));