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:
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.rmsd
program read all the frames from a trajectory, and compute the
root mean square displacement of the 10th atom in the frame.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.#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;
}
#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;
}
#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;
}
#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;
}
#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;
}
#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;
}