# Analyzing a lammps trajectory¶

In this example, a lammps trajectory in dump-text format will be read in, and Steinhardt’s parameters will be calculated.

:

import pyscal as pc
import os
import pyscal.traj_process as ptp
import matplotlib.pyplot as plt
import numpy as np


First, we will use the split_trajectory method from pyscal.traj_process module to help split the trajectory into individual snapshots.

:

trajfile = "traj.light"
files = ptp.split_trajectory(trajfile)


files contain the individual time slices from the trajectory.

:

len(files)

:

10

:

files

:

'traj.light.snap.0.dat'


Now we can make a small function which reads a single configuration and calculates $$q_6$$ values.

:

def calculate_q6(file, format="lammps-dump"):
sys = pc.System()
sys.find_neighbors(method="cutoff", cutoff=0)
sys.calculate_q(6)
q6 = sys.get_qvals(6)
return q6


There are a couple of things of interest in the above function. The find_neighbors method finds the neighbors of the individual atoms. Here, an adaptive method is used, but, one can also use a fixed cutoff or Voronoi tessellation. Also only the unaveraged $$q_6$$ values are calculated above. The averaged ones can be calculate using the averaged=True keyword in both calculate_q and get_qvals method. Now we can simply call the function for each file..

:

q6s = [calculate_q6(file) for file in files]


We can now visualise the calculated values

:

plt.plot(np.hstack(q6s))

:

[<matplotlib.lines.Line2D at 0x7f5e424c0f60>] We will now modify the above function to also find clusters which satisfy particular $$q_6$$ value. But first, for a single file.

:

sys = pc.System()
sys.find_neighbors(method="cutoff", cutoff=0)
sys.calculate_q(6)


Now a clustering algorithm can be applied on top using the cluster_atoms method. cluster_atoms uses a condition as argument which should give a True/False value for each atom. Lets define a condition.

:

def condition(atom):
return atom.get_q(6) > 0.5


The above function returns True for any atom which has a $$q_6$$ value greater than 0.5 and False otherwise. Now we can call the cluster_atoms method.

:

sys.cluster_atoms(condition)

:

16


The method returns 16, which here is the size of the largest cluster of atoms which have $$q_6$$ value of 0.5 or higher. If information about all clusters are required, that can also be accessed.

:

atoms = sys.atoms


atom.cluster gives the number of the cluster that each atom belongs to. If the value is -1, the atom does not belong to any cluster, that is, the clustering condition was not met.

:

clusters = [atom.cluster for atom in atoms if atom.cluster != -1]


Now we can see how many unique clusters are there, and what their sizes are.

:

unique_clusters, counts = np.unique(clusters, return_counts=True)


counts contain all the necessary information. len(counts) will give the number of unique clusters.

:

plt.bar(range(len(counts)), counts)
plt.ylabel("Number of atoms in cluster")
plt.xlabel("Cluster ID")

:

Text(0.5, 0, 'Cluster ID') Now we can finally put all of these together into a single function and run it over our individual time slices.

:

def calculate_q6_cluster(file, cutoff_q6 = 0.5, format="lammps-dump"):
sys = pc.System()
sys.find_neighbors(method="cutoff", cutoff=0)
sys.calculate_q(6)
def _condition(atom):
return atom.get_q(6) > cutoff_q6
sys.cluster_atoms(condition)
atoms = sys.atoms
clusters = [atom.cluster for atom in atoms if atom.cluster != -1]
unique_clusters, counts = np.unique(clusters, return_counts=True)
return counts

:

q6clusters = [calculate_q6_cluster(file) for file in files]


We can plot the number of clusters for each slice

:

plt.plot(range(len(q6clusters)), [len(x) for x in q6clusters], 'o-')
plt.xlabel("Time slice")
plt.ylabel("number of unique clusters")

:

Text(0, 0.5, 'number of unique clusters') We can also plot the biggest cluster size

:

plt.plot(range(len(q6clusters)), [max(x) for x in q6clusters], 'o-')
plt.xlabel("Time slice")
plt.ylabel("Largest cluster size")

:

Text(0, 0.5, 'Largest cluster size') The final thing to do is to remove the split files after use.

:

for file in files:
os.remove(file)


## Using ASE¶

The above example can also done using ASE. The ASE read method needs to be imported.

:

from ase.io import read

:

traj = read("traj.light", format="lammps-dump-text", index=":")


In the above function, index=":" tells ase to read the complete trajectory. The individual slices can now be accessed by indexing.

:

traj

:

Atoms(symbols='H500', pbc=True, cell=[18.21922, 18.22509, 18.36899], momenta=...)


We can use the same functions as above, but by specifying a different file format.

:

q6clusters_ase = [calculate_q6_cluster(x, format="ase") for x in traj]


We will plot and compare with the results from before,

:

plt.plot(range(len(q6clusters_ase)), [max(x) for x in q6clusters_ase], 'o-')
plt.xlabel("Time slice")
plt.ylabel("Largest cluster size")

:

Text(0, 0.5, 'Largest cluster size') As expected, the results are identical for both calculations!