Compute virial from a trajectory file in a post-processing step

compute_virial.sh

#! /usr/bin/env python

from __future__ import division, print_function

import sys

from ase.io import NetCDFTrajectory

from atomistica import TabulatedAlloyEAM, KumagaiScr

###

#calc = TabulatedAlloyEAM(fn='Au-Grochola-JCP05.eam.alloy')
calc = KumagaiScr()

###

traj = NetCDFTrajectory('traj.nc', types_to_numbers=[-1,14], mode='a')

for frame, a in enumerate(traj):
    sys.stdout.write('=== {}/{} ==='.format(frame+1, len(traj)))
    a.set_calculator(calc)
    virial = a.get_stresses()
    a.set_array('virial', virial)
    # Dump virial array to NetCDF file.
    traj.write_arrays(a, frame, ['virial'])