Saving Array Quantities#



  • How can I save per-particle quantities for later analysis?

  • How can I access that data?


  • Show how to log per-particle properties to a GSD file.

  • Explain how to read logged quantities from a GSD file.

  • Mention that OVITO reads these quantities.

Boilerplate code#

import gsd.hoomd
import hoomd
import matplotlib
import numpy

%matplotlib inline'ggplot')
import matplotlib_inline


The render function in the next (hidden) cell will render the system state using fresnel. Find the source in the hoomd-examples repository.

General array quantities#

You logged scalar and array quantities to a HDF5 file in the previous section. HDF5 is the ideal format to log this type of data, as it produces small files that can be read quickly.

Array quantities coupled with a trajectory#

You can also log scalar, array, per-particle (e.g. energy and force), and other quantities in a GSD file along with the trajectory. Then you can use utilize the data in your analysis and visualization workflow, such as using OVITO to color particles by energy or display force vectors. When using OVITO, open the GSD file and all logged quantities will be available in the inspector and for use in the pipeline.

Define the Simulation#

This tutorial executes the Lennard-Jones particle simulation from a previous tutorial. See Introducing Molecular Dyamics for a complete description of this code.

cpu = hoomd.device.CPU()
simulation = hoomd.Simulation(device=cpu, seed=1)

integrator =
cell =
lj =
lj.params[('A', 'A')] = dict(epsilon=1, sigma=1)
lj.r_cut[('A', 'A')] = 2.5
nvt =
simulation.operations.integrator = integrator

Logging per-particle quantities#

MD forces provide a number of loggable quantities including their contribution to the system energy, but also per-particle energy contributions (in energies) and per-particle forces, torques, and virials.

{'energy': 'scalar',
 'energies': 'particle',
 'additional_energy': 'scalar',
 'forces': 'particle',
 'torques': 'particle',
 'virials': 'particle',
 'additional_virial': 'sequence'}

Add the per-particle LJ energies and forces to a logger:

logger = hoomd.logging.Logger()
logger.add(lj, quantities=['energies', 'forces'])

Writing per-particle quantities to a GSD file#

Create the GSD writer to write the simulation trajectory:

gsd_writer = hoomd.write.GSD(filename='trajectory.gsd',

Set the logger attribute and GSD will also store the selected quantities.

gsd_writer.logger = logger

Run the simulation:


Flush the GSD buffer so it is readable in the next code cell. In typical workflows, you will run separate simulation and analysis scripts and the buffered writes will automatically flush when your simulation script exits.


Reading logged data from a GSD file#

Use the gsd package to open the file:

traj ='trajectory.gsd', mode='r')

The log data for a specific frame is stored in the log dictionary for that frame.

dict_keys(['particles/md/pair/LJ/energies', 'particles/md/pair/LJ/forces'])

GSD prepends particles/ to the logged name of per-particle quantities. The quantities are NumPy arrays with N_particles elements. Here are a few slices:

array([-2.48199661, -2.12335126, -4.00387093, -2.63144287, -2.18554869,
       -2.32050608, -2.02647941, -2.44569987, -2.53925901, -2.40447241])
array([[ -2.79783236,  -1.1447838 ,  -2.29953834],
       [  1.34913891,   2.73176148,  -2.13988326],
       [  3.18757211, -12.6700367 ,  -1.80576571],
       [ -4.70877614,   3.40071066,  -8.78080297],
       [  4.94235545,   6.16337821,  -6.97153234],
       [ -3.2043855 ,   0.18160722,  -1.14755586],
       [-10.73489983,   0.28556186,   4.2747685 ],
       [ -2.20080722,  -1.70941385,  16.53612282],
       [  4.01219479,   1.53498589,   3.94355533],
       [  2.1649503 ,   1.18676842,  -6.51449987]])

You can use these arrays as inputs to any computation or plotting tools:

numpy.mean(traj[-1].log['particles/md/pair/LJ/forces'], axis=0)
array([ 5.03069808e-17,  1.64798730e-16, -2.04697370e-16])
fig = matplotlib.figure.Figure(figsize=(5, 3.09))
ax = fig.add_subplot()
ax.hist(traj[-1].log['particles/md/pair/LJ/energies'], 100)
ax.set_xlabel('potential energy')

As with scalar quantities, the array quantities are stored separately in each frame. Use a loop to access a range of frames and compute time-series data or averages.

This is what the system looks like when you color the particles by the values in particles/md/pair/LJ/energies:


In this section, you have logged per-particle quantities to a GSD file during a simulation run and accessed that data with a script. The next section of this tutorial demonstrates how to log particle shape information that OVITO can use.