Add method to write new snapshot that contains a subset of another snapshot's molecules.

[chrisjonesBSU]

This method takes in a "bulk" system snapshot, and a list of arrays of particle indices, and creates a new snapshot that contains only the molecules identified by the list of particles. This might be useful for something like, taking a polymer melt, picking out 2 chains from the melt, and using that as input for another use case (chain entanglement analysis, another MD simulation, etc..). This could be useful for welding/diffusion analysis where we want to iterate over all "left" chains and perform some measurement or analysis over all "right" chains.

Here is an example where the input trajectory1.gsd contains a melt of 50 polymers of 50 monomers each, and a new snapshot is created that contains only chains 0 and 20:

import numpy as np
import gsd.hoomd

from cmeutils.gsd_utils import get_molecule_cluster, trim_snapshot_molecules

residue1 = 0
residue2 = 20
residue_clusters = get_molecule_cluster(gsd_file="trajectory1.gsd")
residue1_indices = np.where(residue_clusters == residue1)[0]
residue2_indices = np.where(residue_clusters == residue2)[0]

with gsd.hoomd.open("trajectory1.gsd", "r")  as traj:
    parent_snap = traj[-1]

new_snap = trim_snapshot_molecules(parent_snapshot=parent_snap, mol_indices=[residue1_indices, residue2_indices])

Here is an example of iterating over an entire trajectory and writing out a new trajectory:

import numpy as np
import gsd.hoomd

from cmeutils.gsd_utils import get_molecule_cluster, trim_snapshot_molecules

residue1 = 0
residue2 = 20
residue_clusters = get_molecule_cluster(gsd_file="trajectory1.gsd")
residue1_indices = np.where(residue_clusters == residue1)[0]
residue2_indices = np.where(residue_clusters == residue2)[0]

with gsd.hoomd.open("trimmed-traj.gsd", "w") as new_traj:
    with gsd.hoomd.open("trajectory1.gsd", "r") as old_traj:
        for snap in old_traj:
            new_snap = trim_snapshot_molecules(parent_snapshot=snap, mol_indices=[residue1_indices, residue2_indices])
            new_traj.append(new_snap)

[codecov]

Codecov Report

Attention: Patch coverage is 1.85185% with 53 lines in your changes are missing coverage. Please review.

Project coverage is 90.29%. Comparing base (3419258) to head (51c015d).

Additional details and impacted files

@@            Coverage Diff             @@
##           master      #85      +/-   ##
==========================================
- Coverage   96.33%   90.29%   -6.04%     
==========================================
  Files           8        8              
  Lines         791      845      +54     
==========================================
+ Hits          762      763       +1     
- Misses         29       82      +53     

Files	Coverage Δ
cmeutils/gsd_utils.py	75.17% <1.85%> (-16.50%)	⬇️

[chrisjonesBSU]

I still need to test this out (e.g. taking one of these snapshots and trying to start a simulation in flower), as well as add a unit test.

[erjank]

It might be worth digging into the GSD documentation or asking on their email list whether there are some slicing or selection operations that accomplish this aim without file i/o. If one can iterate through GSD frames selecting only two chains with existing syntax then we needn't build a whole new thing.