New RDF method with flexible pair exclusions defined by an arbitrary bond depth.

cmeutils/__init__.py

@@ -1,4 +1,10 @@
 from . import gsd_utils, structure
 from .__version__ import __version__
 
-__all__ = ["__version__", "gsd_utils", "structure", "dynamics", "visualize"]
+__all__ = [
+    "__version__",
+    "gsd_utils",
+    "structure",
+    "dynamics",
+    "visualize",
+]

cmeutils/__version__.py

@@ -1 +1 @@
-__version__ = "1.4.0"
+__version__ = "1.5.0"

cmeutils/geometry.py

@@ -1,5 +1,49 @@
 import numpy as np
 from numpy.linalg import svd
+from rowan import vector_vector_rotation
+
+
+def get_quaternions(n_views=20):
+    """Generate quaternions over the specified number of views.
+
+    The first (n_view - 3) views will be the views even distributed on a sphere,
+    while the last three views will be the face-on, edge-on, and corner-on
+    views, respectively.
+
+    These quaternions are useful as input to `view_orientation` kwarg in
+    `freud.diffraction.Diffractometer.compute`.
+
+    Parameters
+    ----------
+    n_views : int, default 20
+        The number of views to compute.
+
+    Returns
+    -------
+    list of numpy.ndarray
+        Quaternions as (4,) arrays.
+    """
+    if n_views <= 3 or not isinstance(n_views, int):
+        raise ValueError("Please set n_views to an integer greater than 3.")
+    # Calculate points for even distribution on a sphere
+    ga = np.pi * (3 - 5**0.5)
+    theta = ga * np.arange(n_views - 3)
+    z = np.linspace(1 - 1 / (n_views - 3), 1 / (n_views - 3), n_views - 3)
+    radius = np.sqrt(1 - z * z)
+    points = np.zeros((n_views, 3))
+    points[:-3, 0] = radius * np.cos(theta)
+    points[:-3, 1] = radius * np.sin(theta)
+    points[:-3, 2] = z
+
+    # face on
+    points[-3] = np.array([0, 0, 1])
+    # edge on
+    points[-2] = np.array([0, 1, 1])
+    # corner on
+    points[-1] = np.array([1, 1, 1])
+
+    unit_z = np.array([0, 0, 1])
+    return [vector_vector_rotation(i, unit_z) for i in points]
 
 
 def get_backbone_vector(coordinates):

cmeutils/gsd_utils.py

@@ -8,6 +8,24 @@
 from boltons.setutils import IndexedSet
 
 
+def snapshot_to_graph(snap):
+    """Convert a HOOMD snapshot into a networkx bond graph.
+
+    Parameters
+    ----------
+    snap : gsd.hoomd.Frame, required
+        Snapshot (frame) containing particle positions and bond groups.
+    """
+    graph = nx.Graph()
+    positions = snap.particles.position
+    type_ids = snap.particles.typeid
+    types = snap.particles.types
+    for i in range(snap.particles.N):
+        graph.add_node(i, type=types[type_ids[i]], pos=positions[i])
+    graph.add_edges_from(snap.bonds.group)
+    return graph
+
+
 def frame_to_freud_system(frame, ref_length=None):
     """Creates a freud system given a gsd.hoomd.Frame.
 
@@ -335,6 +353,50 @@ def identify_snapshot_connections(snapshot):
     return snapshot
 
 
+def get_centers(gsdfile, new_gsdfile):
+    """Create a gsd file containing the molecule centers from an existing gsd
+     file.
+
+
+    This function calculates the centers of a trajectory given a GSD file
+    and stores them into a new GSD file just for centers. By default it will
+    calculate the centers of an entire trajectory.
+
+    Parameters
+    ----------
+    gsdfile : str
+        Filename of the GSD trajectory.
+    new_gsdfile : str
+        Filename of new GSD for centers.
+    """
+    with (
+        gsd.hoomd.open(new_gsdfile, "w") as new_traj,
+        gsd.hoomd.open(gsdfile, "r") as traj,
+    ):
+        snap = traj[0]
+        cluster_idx = get_molecule_cluster(snap=snap)
+        for snap in traj:
+            new_snap = gsd.hoomd.Frame()
+            new_snap.configuration.box = snap.configuration.box
+            f_box = freud.box.Box.from_box(snap.configuration.box)
+            # Use the freud box to unwrap the particle positions
+            unwrapped_positions = f_box.unwrap(
+                snap.particles.position, snap.particles.image
+            )
+            uw_centers = []
+            for i in range(max(cluster_idx) + 1):
+                cluster_uw_pos = unwrapped_positions[np.where(cluster_idx == i)]
+                uw_centers.append(np.mean(cluster_uw_pos, axis=0))
+            uw_centers = np.stack(uw_centers)
+            new_snap.particles.position = f_box.wrap(uw_centers)
+            new_snap.particles.N = len(uw_centers)
+            new_snap.particles.types = ["A"]
+            new_snap.particles.image = f_box.get_images(uw_centers)
+            new_snap.particles.typeid = np.zeros(len(uw_centers))
+            new_snap.validate()
+            new_traj.append(new_snap)
+
+
 def _fill_connection_info(snapshot, connections, type_):
     p_types = snapshot.particles.types
     p_typeid = snapshot.particles.typeid