Usage examples

Read a single frame (indexes.py)

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from chemfiles import Trajectory

trajectory = Trajectory("filename.xyz")
frame = trajectory.read()
positions = frame.positions()

indexes = []
for i in range(frame.natoms()):
    # positions is a numpy ndarray
    if positions[i, 0] < 5:
        indexes.append(i)

print("Atoms with x < 5:")
for i in indexes:
    print("  - {}".format(i))

Read Multiple frames (rmsd.py)

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from chemfiles import Trajectory
import math

trajectory = Trajectory("filename.nc")
distances = []
# Accumulate the distances to the origin of the 10th atom throughtout the
# trajectory
for frame in trajectory:
    # Position of the 10th atom
    position = frame.positions()[9, :]
    distance = math.sqrt(position.dot(position))
    distances.append(distance)

mean = sum(distances) / len(distances)
rmsd = 0.0
for dist in distances:
    rmsd += (mean - dist) * (mean - dist)

rmsd /= len(distances)
rmsd = math.sqrt(rmsd)

print("Root-mean square displacement is: {}".format(rmsd))

Write frames (convert.py)

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from chemfiles import Trajectory, UnitCell, Atom, Topology

trajectory = Trajectory("water.xyz")

# Orthorombic UnitCell with lengths of 20, 15 and 35 A
trajectory.set_cell(UnitCell(20, 15, 35))

# Create topology with one water molecule
topology = Topology()
topology.add_atom(Atom("O"))
topology.add_atom(Atom("H"))
topology.add_atom(Atom("H"))

topology.add_bond(0, 1)
topology.add_bond(0, 2)
trajectory.set_topology(topology)

# Write an output file using PDB format
with Trajectory("water.pdb", "w") as output:
    for frame in trajectory:
        output.write(frame)

trajectory.close()

Use selections (select.py)

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from chemfiles import Trajectory, Selection

# Read a frame from a given file
frame = Trajectory("filename.xyz").read()


# Create a selection for all atoms with "Zn" name
selection = Selection("name Zn")
# Get the list of matching atoms from the frame
zincs = selection.evaluate(frame)

print("We have {} zinc in the frame".format(len(zincs)))
for i in zincs:
    print("{} is a zinc".format(i))


# Create a selection for multiple atoms
selection = Selection("angles: name(#1) H and name(#2) O and name(#3) H")
# Get the list of matching atoms in the frame
waters = selection.evaluate(frame)

print("We have {} water molecules".format(len(waters)))
for (i, j, k) in waters:
    print("{} - {} - {} is a water".format(i, j, k))