Usage examples

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)

This page contains the following examples, written in all the languages supported by chemfiles:

  • The indexes program read the first frame of a trajectory and fill a vector called indexes with the indexes of all the atom obeying the \(x < 5\) relation.
  • The rmsd program read all the frames from a trajectory, and compute the root mean square displacement of the 10th atom in the frame.
  • The convert program read all the frames from a trajectory, set the unit cell and the topology of the frame, and then write it back to another file. The input file contain here the trajectory of a single water molecule, but without any topological or cell information. So the program add these informations to a new trajectory file.

Usage from C++

Read a single frame (indexes.cpp)

#include <iostream>
#include <vector>

#include "chemfiles.hpp"

int main() {
    chemfiles::Trajectory file("filename.xyz");

    auto frame = file.read();
    auto positions = frame.positions();
    std::vector<size_t> indexes;

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

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

    return 0;
}

Read Multiple frames (rmsd.cpp)

#include <iostream>
#include <vector>
#include <numeric>
#include <cmath>

#include "chemfiles.hpp"

int main() {
    chemfiles::Trajectory file("filename.nc");
    std::vector<double> distances;

    // Accumulate the distances to the origin of the 10th atom throughtout the
    // trajectory
    chemfiles::Frame frame;
    while (!file.done()) {
        frame = file.read();
        // Position of the 10th atom
        auto pos = frame.positions()[9];
        double distance = sqrt(pos[0]*pos[0] + pos[1]*pos[1] + pos[2]*pos[2]);
        distances.push_back(distance);
    }

    double mean = std::accumulate(std::begin(distances), std::end(distances), 0.0) / static_cast<double>(distances.size());
    double rmsd = 0.0;
    for (auto dist : distances) {
        rmsd += (mean - dist)*(mean - dist);
    }
    rmsd /= static_cast<double>(distances.size());
    rmsd = sqrt(rmsd);

    std::cout << "Root-mean square displacement is: " << rmsd << std::endl;
    return 0;
}

Write frames (convert.cpp)

#include "chemfiles.hpp"

int main() {
    // Open the input file for reading
    chemfiles::Trajectory input("water.xyz");

    // Set the unit cell to an orthorombic cell with lengths of 20, 15 and 35 A
    input.set_cell(chemfiles::UnitCell(20, 15, 35));

    // Create a water molecule topology
    chemfiles::Topology water;
    water.append(chemfiles::Atom("O"));
    water.append(chemfiles::Atom("H"));
    water.append(chemfiles::Atom("H"));
    water.add_bond(0, 1);
    water.add_bond(0, 2);
    input.set_topology(water);

    // Write to the output file using PDB format
    chemfiles::Trajectory output("water.pdb", 'w');
    while (!input.done()) {
        // The UnitCell and the Topology are automatically set when reading a
        // frame to the specified cell and topology.
        auto frame = input.read();
        output.write(frame);
    }

    return 0;
}

Usage from C

Read a single frame (indexes.c)

#include <stdio.h>
#include <stdlib.h>

#include "chemfiles.h"

int main(void) {
    CHFL_TRAJECTORY* file = chfl_trajectory_open("tests/files/xyz/helium.xyz", 'r');
    CHFL_FRAME* frame = chfl_frame();
    unsigned* indexes = NULL;

    if (chfl_trajectory_read(file, frame) != CHFL_SUCCESS) {/*Handle error*/}

    uint64_t natoms = 0;
    chfl_vector_t* positions = NULL;
    chfl_frame_positions(frame, &positions, &natoms);
    indexes = malloc((size_t)natoms * sizeof(unsigned));
    if (indexes == NULL) {/*Handle error*/}

    for (unsigned i=0; i<natoms; i++) {
        indexes[i] = (unsigned)-1;
    }

    unsigned last_index = 0;
    for (unsigned i=0; i<natoms; i++) {
        if (positions[i][0] < 5) {
            indexes[last_index] = i;
            last_index++;
        }
    }

    printf("Atoms with x < 5:\n");
    unsigned i = 0;
    while(indexes[i] != (unsigned)-1 && i < natoms) {
        printf("  - %d\n", indexes[i]);
        i++;
    }
    printf("Number of atoms: %d\n", i);

    chfl_trajectory_close(file);
    chfl_frame_free(frame);
    free(indexes);
    return EXIT_SUCCESS;
}

Read Multiple frames (rmsd.c)

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "chemfiles.h"

int main(void) {
    CHFL_TRAJECTORY* file = chfl_trajectory_open("filename.nc", 'r');
    uint64_t nsteps = 0;
    chfl_trajectory_nsteps(file, &nsteps);

    double* distances = malloc((size_t)nsteps * sizeof(double));
    if (distances == NULL) {/*Handle error*/}

    CHFL_FRAME* frame = chfl_frame();
    // Accumulate the distances to the origin of the 10th atom throughtout the
    // trajectory
    for (uint64_t i=0; i<nsteps; i++) {
        if(!chfl_trajectory_read(file, frame)) {/*Handle error*/}

        uint64_t natoms = 0;
        chfl_vector_t* positions = NULL;

        // Get a pointer to the positions in `positions`. The array `positions`
        // contains natoms entries.
        chfl_frame_positions(frame, &positions, &natoms);
        if (natoms < 10) {
            printf("Not enough atoms in the frame\n");
            continue;
        }
        // Get the distance to the origin of the 10th atom
        double distance = sqrt(positions[9][0] * positions[9][0] +
                               positions[9][1] * positions[9][1] +
                               positions[9][2]*positions[9][2]);
        distances[i] = distance;
    }

    double mean = 0;
    for (uint64_t i=0; i<nsteps; i++) {
        mean += distances[i];
    }
    mean /= (double)nsteps;

    double rmsd = 0.0;
    for (uint64_t i=0; i<nsteps; i++) {
        rmsd += (mean - distances[i])*(mean - distances[i]);
    }
    rmsd /= (double)nsteps;
    rmsd = sqrt(rmsd);

    printf("Root-mean square displacement is: %f", rmsd);

    chfl_trajectory_close(file);
    chfl_frame_free(frame);
    free(distances);
    return 0;
}

Write frames (convert.c)

#include <stdio.h>
#include <stdlib.h>

#include "chemfiles.h"

int main(void) {
    CHFL_TRAJECTORY* input = chfl_trajectory_open("water.xyz", 'r');

    // Set the unit cell to an orthorombic cell with lengths of 20, 15 and 35 A
    CHFL_CELL* cell = chfl_cell((chfl_vector_t){20, 15, 35});
    chfl_trajectory_set_cell(input, cell);
    chfl_cell_free(cell);

    // Create a water molecule topology
    CHFL_ATOM* O = chfl_atom("O");
    CHFL_ATOM* H = chfl_atom("H");
    CHFL_TOPOLOGY* water = chfl_topology();
    chfl_topology_add_atom(water, O);
    chfl_topology_add_atom(water, H);
    chfl_topology_add_atom(water, H);
    chfl_topology_add_bond(water, 0, 1);
    chfl_topology_add_bond(water, 0, 2);

    chfl_trajectory_set_topology(input, water);

    chfl_topology_free(water);
    chfl_atom_free(O);
    chfl_atom_free(H);

    uint64_t nsteps = 0;
    chfl_trajectory_nsteps(input, &nsteps);

    CHFL_TRAJECTORY* output = chfl_trajectory_open("water.pdb", 'w');
    CHFL_FRAME* frame = chfl_frame();
    for (uint64_t i=0; i<nsteps; i++) {
        // The unit cell and the topology are automatically set when reading a
        // frame.
        chfl_trajectory_read(input, frame);
        chfl_trajectory_write(output, frame);
    }

    chfl_trajectory_close(input);
    chfl_trajectory_close(output);
    chfl_frame_free(frame);
    return 0;
}