Extraction of cluster geometries from a trajectory

from exa.util import isotopes
import exatomic
from exatomic.core.two import compute_atom_two_out_of_core    # If we need to compute atom_two out of core (low RAM)
from exatomic.algorithms.neighbors import periodic_nearest_neighbors_by_atom    # Only valid for simple cubic periodic cells
from exatomic.base import resource # Helper function to find the path to an example file

xyzpath = resource("rh2.traj.xyz")

xyz = exatomic.XYZ(xyzpath)    # Editor that can parse the XYZ file and build a universe!

u = xyz.to_universe()    # Parse the XYZ data to a universe

u.info()    # Current tables

	size	type
object
METADATA	48	-
WIDGET	0	-
atom	985706	exatomic.core.atom.Atom
frame	1616	exatomic.core.frame.Frame

# Add the unit cell dimensions to the frame table of the universe
a = 37.08946302
for i, q in enumerate(("x", "y", "z")):
    for j, r in enumerate(("i", "j", "k")):
        if i == j:
            u.frame[q+r] = a
        else:
            u.frame[q+r] = 0.0
    u.frame["o"+q] = 0.0
u.frame['periodic'] = True

len(u)    # Number of steps in this trajectory (small for example purposes)

101

isotopes.Rh.radius    # Default Rh (covalent) radius

2.3621574999999995

To extract clusters from the trajectory, we don’t need to pre-compute atom_two

u.atom.head()

	symbol	x	y	z	frame
atom
0	Rh	19.3251	20.99490	12.9253	0
1	Rh	20.8605	19.40470	17.0475	0
2	Cl	32.7903	9.48909	27.4356	0
3	Cl	27.2017	10.59920	26.4768	0
4	Cl	39.1224	25.75620	13.6518	0

%%time
dct = periodic_nearest_neighbors_by_atom(u,       # Universe with all frames from which we want to extract clusters
                                         "Rh",    # Source atom from which we will search for neighbors
                                         a,       # Cubic cell dimension
                                         # Cluster sizes we want
                                         [0, 1, 2, 3, 4, 8, 12, 16],
                                         # Additional arguments to be passed to the atomic two body calculation
                                         # Note that it is critical to increase dmax (in compute_atom_two)!!!
                                         Rh=2.6, dmax=a/2)

CPU times: user 2min 20s, sys: 2min 35s, total: 4min 56s
Wall time: 5min 20s

dct.keys()    # 'dct' is a dictionary containing our clustered universes as well as the sorted neighbor list (per frame)

dict_keys(['nearest', 0, 1, 2, 3, 4, 8, 12, 16])

dct['nearest'].head()    # Table of sorted indexes of nearest molecules to "Rh" in each frame

	frame	idx	two	atom	molecule
0	0	18	421055	5186	890
1	0	26	420722	4930	879
2	0	30	420997	5083	859
3	0	34	420801	5090	857
4	0	35	420834	5155	881

exatomic.UniverseWidget(dct[0])    # Just the analyte

exatomic.UniverseWidget(dct[16])   # View the universe with 16 nearest molecules

If we had wanted to view the entire universe…

%time u.compute_atom_two(Rh=2.6)    # Compute atom two body data with slightly larger Rh radius
# OR
#%time compute_atom_two_out_of_core("atom_two.hdf", u, a, Rh=2.6)      # Writes the two body data to "atom_two.hdf" with keys "frame_{index}/atom_two"

CPU times: user 3.17 s, sys: 281 ms, total: 3.45 s
Wall time: 3.51 s

u.memory_usage().sum()/1024**2    # RAM usage (MB)

10.53761100769043

exatomic.UniverseWidget(u)