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

<narupatools.frame.hdf5._traj.HDF5Trajectory at 0x7f56902759a0>

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.07722139, 3.43267345, 3.7069149 ])

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

1262.3096923828125

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

1.2569999694824219

# Maximum kinetic energy achieved
max(trajectory.kinetic_energies)

2707.7802734375

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())

[<narupatools.frame.hdf5._traj.HDF5Interaction at 0x7f565a128220>,
 <narupatools.frame.hdf5._traj.HDF5Interaction at 0x7f565a1285e0>,
 <narupatools.frame.hdf5._traj.HDF5Interaction at 0x7f565a1288b0>,
 <narupatools.frame.hdf5._traj.HDF5Interaction at 0x7f565a1289d0>,
 <narupatools.frame.hdf5._traj.HDF5Interaction at 0x7f565a128a90>,
 <narupatools.frame.hdf5._traj.HDF5Interaction at 0x7f565a128b50>,
 <narupatools.frame.hdf5._traj.HDF5Interaction at 0x7f565a128bb0>]

We can access individual interactions and look at their properties

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

<narupatools.frame.hdf5._traj.HDF5Interaction at 0x7f565a128d30>

# Type used for the interaction
interaction.interaction_type

'spring'

# Duration of the interaction
interaction.duration

0.012000084

# Particle indices affected by this interaction
interaction.particle_indices

array([60])

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

array([5000., 5000., 5000., 5000., 5000., 5000., 5000., 5000., 5000.,
       5000., 5000., 5000., 5000.])

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.0609578669076

interaction.potential_energies

array([  1.47560036,   1.71617949,   0.96935296,   0.82653397,
         3.84148312,   9.00822926,  36.86339951,  52.28302002,
        85.17850494, 123.67550659, 147.31652832, 164.91673279,
       211.69317627])

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 import show_trajectory

show_trajectory(trajectory)