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 traj("filename.xyz");
    chemfiles::Frame frame;

    traj >> frame;
    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 traj("filename.nc");
    chemfiles::Frame frame;

    std::vector<double> distances;

    // Accumulate the distances to the origin of the 10th atom throughtout the
    // trajectory
    while (!traj.done()) {
        traj >> frame;
        // 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() {
    chemfiles::Trajectory input("water.xyz");
    chemfiles::Frame frame{};
    chemfiles::Topology water_topology{};
    // Orthorombic UnitCell with lengths of 20, 15 and 35 A
    chemfiles::UnitCell cell(20, 15, 35);

    // Create Atoms
    chemfiles::Atom O("O");
    chemfiles::Atom H("H");

    // Fill the topology with one water molecule
    water_topology.append(O);
    water_topology.append(H);
    water_topology.append(H);
    water_topology.add_bond(0, 1);
    water_topology.add_bond(0, 2);

    chemfiles::Trajectory output("water.pdb", "w");

    while (!input.done()) {
        input >> frame;
        // Set the frame cell and topology
        frame.set_cell(cell);
        frame.set_topology(water_topology);
        // Write the frame to the output file, using PDB format
        output << frame;
    }

    return 0;
}

Usage from C

Read a single frame (indexes.c)

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

#include "chemfiles.h"

int main() {
    CHFL_TRAJECTORY* traj = chfl_trajectory_open("tests/files/xyz/helium.xyz", "r");
    CHFL_FRAME* frame = chfl_frame(0);

    if (traj == NULL)
        goto error;

    if (chfl_trajectory_read(traj, frame))
        goto error;

    size_t natoms = 0;
    float (*positions)[3] = NULL;
    chfl_frame_positions(frame, &positions, &natoms);
    unsigned* indexes = (unsigned*)malloc(natoms*sizeof(unsigned));
    if (indexes == NULL)
        goto error;

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

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

    printf("Atoms with x < 5:\n");
    int 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(traj);
    chfl_frame_free(frame);
    free(positions);
    return 0;

error:
    printf("Error, cleaning up …\n");
    chfl_trajectory_close(traj);
    chfl_frame_free(frame);
    free(positions);
    return 1;
}

Read Multiple frames (rmsd.c)

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

#include "chemfiles.h"

int main() {
    CHFL_TRAJECTORY* traj = chfl_trajectory_open("filename.nc", "r");
    CHFL_FRAME* frame = chfl_frame(0);
    float (*positions)[3] = NULL;
    double* distances = NULL;

    if(traj == NULL || frame == NULL)
        goto error;

    size_t nsteps = 0;
    chfl_trajectory_nsteps(traj, &nsteps);

    distances = (double*)malloc(sizeof(double)*nsteps);
    if (distances == NULL)
        goto error;

    // Accumulate the distances to the origin of the 10th atom throughtout the
    // trajectory
    for (size_t i=0; i<nsteps; i++) {
        if(!chfl_trajectory_read(traj, frame))
            goto error;

        size_t natoms = 0;
        // Position of the 10th atom
        chfl_frame_positions(frame, &positions, &natoms);
        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 (size_t i=0; i<nsteps; i++) {
        mean += distances[i];
    }
    mean /= nsteps;

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

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

    // Free the memory
    chfl_trajectory_close(traj);
    chfl_frame_free(frame);
    free(distances);
    free(positions);
    return 0;

error:
    printf("Error, cleaning up …\n");
    chfl_trajectory_close(traj);
    chfl_frame_free(frame);
    free(distances);
    free(positions);
    return 1;
}

Write frames (convert.c)

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

#include "chemfiles.h"

int main() {
    CHFL_TRAJECTORY* input = chfl_trajectory_open("water.xyz", "r");
    CHFL_TRAJECTORY* output = chfl_trajectory_open("water.pdb", "w");
    CHFL_FRAME* frame = chfl_frame(0);
    CHFL_TOPOLOGY* water_topology = chfl_topology();
    // Orthorombic UnitCell with lengths of 20, 15 and 35 A
    CHFL_CELL* cell = chfl_cell(20, 15, 35);

    // Create Atoms
    CHFL_ATOM* O = chfl_atom("O");
    CHFL_ATOM* H = chfl_atom("H");

    if (input == NULL || frame == NULL || water_topology == NULL ||
        cell == NULL || O == NULL || H == NULL || output == NULL)
            goto error;

    // Fill the topology with one water molecule
    chfl_topology_append(water_topology, O);
    chfl_topology_append(water_topology, H);
    chfl_topology_append(water_topology, H);
    chfl_topology_add_bond(water_topology, 0, 1);
    chfl_topology_add_bond(water_topology, 0, 2);

    size_t nsteps;
    chfl_trajectory_nsteps(input, &nsteps);

    for (size_t i=0; i<nsteps; i++) {
        chfl_trajectory_read(input, frame);
        // Set the frame cell and topology
        chfl_frame_set_cell(frame, cell);
        chfl_frame_set_topology(frame, water_topology);
        // Write the frame to the output file, using PDB format
        chfl_trajectory_write(output, frame);
    }

    chfl_trajectory_close(input);
    chfl_trajectory_close(output);
    chfl_frame_free(frame);
    chfl_cell_free(cell);
    chfl_topology_free(water_topology);
    chfl_atom_free(O);
    chfl_atom_free(H);
    return 0;

error:
    printf("Error, cleaning up …\n");
    chfl_trajectory_close(input);
    chfl_trajectory_close(output);
    chfl_frame_free(frame);
    chfl_cell_free(cell);
    chfl_topology_free(water_topology);
    chfl_atom_free(O);
    chfl_atom_free(H);
    return 1;
}