Frame

class 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

explicit 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());

inline 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);

inline 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

Throws:

Error – if the topology size does not match the size of this frame

inline 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);

inline 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);

inline 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);

inline span<Vector3D> positions()

Get the positions (in Angstroms) 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));

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

Get the positions (in Angstroms) of the atoms 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());

inline optional<span<Vector3D>> velocities()

Get an velocities (in Angstroms/ps) 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));

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

Get an velocities (in Angstroms/ps) 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));

Throws:

chemfiles::OutOfBounds – if i is bigger than the number of atoms in this frame

inline 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);

inline 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.

// 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);

Throws:

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

inline 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);

inline 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

Throws:

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

inline 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

Throws:

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

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

inline 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

Throws:

OutOfBounds – if atom_i or atom_j are greater than size()

inline 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

Throws:

OutOfBounds – if index is greater than size()

inline 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

Throws:

OutOfBounds – if index is greater than size()

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

Throws:

chemfiles::OutOfBounds – if i or j are bigger than the number of atoms in this frame

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

Throws:

chemfiles::OutOfBounds – if i, j or k are bigger than the number of atoms in this frame

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

Throws:

chemfiles::OutOfBounds – if i, j, k or m are bigger than the number of atoms in this frame

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

Throws:

chemfiles::OutOfBounds – if i, j, k or m are bigger than the number of atoms in this frame

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

inline 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);

inline 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>
inline 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"));