Distinction of solid liquid atoms and clustering
In this example, we will take one snapshot from a molecular dynamics simulation which has a solid cluster in liquid. The task is to identify solid atoms and cluster them. More details about the method can be found here.
The first step is, of course, importing all the necessary module. For visualisation, we will use Ovito.
The above image shows a visualisation of the system using Ovito. Importing modules,
[1]:
import pyscal.core as pc
Now we will set up a System with this input file, and calculate neighbors. Here we will use a cutoff method to find neighbors. More details about finding neighbors can be found here.
[2]:
sys = pc.System()
sys.read_inputfile('cluster.dump')
sys.find_neighbors(method='cutoff', cutoff=3.63)
Once we compute the neighbors, the next step is to find solid atoms. This can be done using System.find_solids
method. There are few parameters that can be set, which can be found in detail here.
[3]:
sys.find_solids(bonds=6, threshold=0.5, avgthreshold=0.6, cluster=False)
The above statement found all the solid atoms. Solid atoms can be identified by the value of the solid
attribute. For that we first get the atom objects and select those with solid
value as True.
[4]:
atoms = sys.atoms
solids = [atom for atom in atoms if atom.solid]
len(solids)
[4]:
203
There are 202 solid atoms in the system. In order to visualise in Ovito, we need to first write it out to a trajectory file. This can be done with the help of to_file
method of System. This method can help to save any attribute of the atom or ant Steinhardt parameter value.
[6]:
sys.to_file('sys.solid.dat', customkeys = ['solid'])
We can now visualise this file in Ovito. After opening the file in Ovito, the modifier compute property can be selected. The Output property
should be selection
and in the expression field, solid==0
can be selected to select all the non solid atoms. Applying a modifier delete selected particles can be applied to delete all the non
solid particles. The system after removing all the liquid atoms is shown below.
Clustering algorithm
You can see that there is a cluster of atom. The clustering functions that pyscal offers helps in this regard. If you used find_clusters
with cluster=True
, the clustering is carried out. Since we did used cluster=False
above, we will rerun the function
[7]:
sys.find_solids(bonds=6, threshold=0.5, avgthreshold=0.6, cluster=True)
[7]:
176
You can see that the above function call returned the number of atoms belonging to the largest cluster as an output. In order to extract atoms that belong to the largest cluster, we can use the largest_cluster
attribute of the atom.
[8]:
atoms = sys.atoms
largest_cluster = [atom for atom in atoms if atom.largest_cluster]
len(largest_cluster)
[8]:
176
The value matches that given by the function. Once again we will save this information to a file and visualise it in Ovito.
[9]:
sys.to_file('sys.cluster.dat', customkeys = ['solid', 'largest_cluster'])
The system visualised in Ovito following similar steps as above is shown below.
It is clear from the image that the largest cluster of solid atoms was successfully identified. Clustering can be done over any property. The following example with the same system will illustrate this.
Clustering based on a custom property
In pyscal, clustering can be done based on any property. The following example illustrates this. To find the clusters based on a custom property, the System.clusters_atoms
method has to be used. The simulation box shown above has the centre roughly at (25, 25, 25). For the custom clustering, we will cluster all atoms within a distance of 10 from the the rough centre of the box at (25, 25, 25). Let us define a function that checks the above condition.
[10]:
def check_distance(atom):
#get position of atom
pos = atom.pos
#calculate distance from (25, 25, 25)
dist = ((pos[0]-25)**2 + (pos[1]-25)**2 + (pos[2]-25)**2)**0.5
#check if dist < 10
return (dist <= 10)
The above function would return True or False depending on a condition and takes the Atom as an argument. These are the two important conditions to be satisfied. Now we can pass this function to cluster. First, set up the system and find the neighbors.
[11]:
sys = pc.System()
sys.read_inputfile('cluster.dump')
sys.find_neighbors(method='cutoff', cutoff=3.63)
Now cluster
[12]:
sys.cluster_atoms(check_distance)
[12]:
242
There are 242 atoms in the cluster! Once again we can check this, save to a file and visualise in ovito.
[13]:
atoms = sys.atoms
largest_cluster = [atom for atom in atoms if atom.largest_cluster]
len(largest_cluster)
[13]:
242
[14]:
sys.to_file('sys.dist.dat', customkeys = ['solid', 'largest_cluster'])
This example illustrates that any property can be used to cluster the atoms!