Reading HDF5 Trajectories

Reading narupatools HDF5 trajectory files is simple.

from narupatools.frame.hdf5 import HDF5Trajectory

# Read in the trajectory with the given filename
trajectory = HDF5Trajectory.load_file("nanotube.h5")
trajectory

Warning: importing 'simtk.openmm' is deprecated.  Import 'openmm' instead.

<HDF5Trajectory 1258 frame(s), 1.257 ps long, 7 interactions(s), 65 atom(s), 1 residues(s), 1 chain(s), 88 bond(s)>

Accessing Trajectory Data

You can access positions, velocities, forces, times, kinetic energies and potential energies from the trajectory

# Position of atom 55 at the end of 43rd frame
trajectory.positions[43][55]

array([4.0772214, 3.4326735, 3.706915 ], dtype=float32)

# Potential energy at the end of the 123rd frame
trajectory.potential_energies[123]

1262.3097

# Duration of the simulation
trajectory.times[-1]

1.257

# Maximum kinetic energy achieved
max(trajectory.kinetic_energies)

2707.7803

Exploring Interactions

Interactions can be viewed through the interactions attribute

trajectory.interactions

<InteractionsView 7 interaction(s)>

This is a dictionary of interaction IDs to interactions

list(trajectory.interactions.keys())

['interaction.166f117c-f73b-4a66-83f6-d1d233158305',
 'interaction.1d63599b-1174-4fbd-b9a2-bcd9af8205b1',
 'interaction.63ec07ab-5b3e-45a1-aa10-f6f7f85decf0',
 'interaction.75546256-6793-4de2-a333-de80da48259b',
 'interaction.a61a910a-6f2c-41eb-8efd-7a80844cc75e',
 'interaction.e4bbb3de-70cd-42b3-9fa8-e0198c1f4b89',
 'interaction.fa99c0a3-a0c6-41b8-a9b1-5ec9aa35764d']

list(trajectory.interactions.values())

[<InteractionView>,
 <InteractionView>,
 <InteractionView>,
 <InteractionView>,
 <InteractionView>,
 <InteractionView>,
 <InteractionView>]

We can access individual interactions and look at their properties

interaction = list(trajectory.interactions.values())[0]
interaction

<InteractionView>

# Type used for the interaction
interaction.type

'spring'

# Duration of the interaction
interaction.duration

0.012000084

# Particle indices affected by this interaction
interaction.indices

array([60])

# Interaction scales for all the frames it was active
interaction.interaction_scales

array([5000., 5000., 5000., 5000., 5000., 5000., 5000., 5000., 5000.,
       5000., 5000., 5000., 5000.], dtype=float32)

interaction.frame_indices

array([1045, 1046, 1047, 1048, 1049, 1050, 1051, 1052, 1053, 1054, 1055,
       1056, 1057])

# Range of frame's this interaction was active
interaction.frame_range

range(1045, 1058)

# Calculate work done by this interaction
interaction.calculate_work() / interaction.duration

-750.0609564988476

interaction.potential_energies

array([  1.4756004 ,   1.7161795 ,   0.96935296,   0.826534  ,
         3.841483  ,   9.008229  ,  36.8634    ,  52.28302   ,
        85.178505  , 123.67551   , 147.31653   , 164.91673   ,
       211.69318   ], dtype=float32)

trajectory.interactions.potential_energies

array([0., 0., 0., ..., 0., 0., 0.])

Plotting Data

The following shows how to use plotly to plot the various energies and work associated with your trajectory.

It uses the calculate_cumulative_work which gives the cumulative work for each frame of the trajectory.

import plotly.graph_objects as go

# Create the figure
fig = go.Figure()

# Add x and y labels
fig.update_yaxes(title_text="Energy (kJ/mol)")
fig.update_xaxes(title_text="Time (ps)")

# Add a kinetic energy line
fig.add_trace(go.Scatter(x=trajectory.times, y=trajectory.kinetic_energies,
                    mode='lines',
                    name='Kinetic Energy'))

# Add a potential energy line, removing the potential energy from interactions
fig.add_trace(go.Scatter(x=trajectory.times, y=trajectory.potential_energies - trajectory.interactions.potential_energies,
                    mode='lines',
                    name='Potential Energy'))

# Add a cumulative work line
fig.add_trace(go.Scatter(x=trajectory.times, y=trajectory.interactions.calculate_cumulative_work(),
                    mode='lines',
                    name='Cumulative Work'))

# Add a interaction energy line
fig.add_trace(go.Scatter(x=trajectory.times, y=trajectory.interactions.potential_energies,
                    mode='lines',
                    name='Interaction Energies'))

# Highlight regions in which interactions were applied
for interaction in trajectory.interactions.values():
    fig.add_vrect(
        x0=interaction.start_time, x1=interaction.end_time,
        fillcolor="LightSalmon", opacity=0.5,
        layer="below", line_width=0,
    )

# Show the figure
fig.show()

NGLView

The trajectory can also be visualised in nglview

from narupatools.nglview.show import show_narupatools_traj
show_narupatools_traj(trajectory)