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.

Array quantities#

You logged a 6-element array quantity in the previous section, the pressure tensor. You can also log larger arrays and per-particle quantities (e.g. energy and force) 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()
sim = 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 =
sim.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#

In the previous section, you used the Null filter to produce GSD file with only logged data. In this section, include all particles so that you can associate logged per-particle quantities with the particle properties.

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

Run the simulation:


As discussed in the previous section, flush the write buffer (this is not necessary in typical workflows).


Reading logged data from a GSD file#

Use the gsd package to open the file:

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

The previous section demonstrated how read_log provides time-series data for plotting and analysis. This section shows how 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.56391048, -1.88822176, -1.41906894, -1.62426695, -2.26188433,
       -1.64826228, -1.86217053, -0.9897792 , -2.7154282 , -1.88687963])
array([[ -2.13808166, -19.36630853,  -8.67495909],
       [ -3.28842103,   7.45784515, -12.86546371],
       [  1.46057136,   0.0920419 ,   0.60696098],
       [  3.43611667,  -2.53969686,   1.88552897],
       [  2.66370616,   7.12540725,   8.06113434],
       [ -3.0789833 ,  -6.1846985 ,  -2.16331169],
       [ -7.16193864,   0.16041416,   3.28611542],
       [-21.36195463,   1.10226678,  26.49873349],
       [ -9.56284366,   0.84839853,   0.09193915],
       [  3.31583343,  -1.10904533,  -1.44974825]])

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([ 1.63064007e-16, -9.71445147e-17, -7.63278329e-17])
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.