Calculating bond orientational order parameters#
This example illustrates the calculation of bond orientational order parameters. Bond order parameters, \(q_l\) and their averaged versions, \(\bar{q}_l\) are widely used to identify atoms belong to different crystal structures. In this example, we will consider bcc, fcc, and hcp, and calculate the \(q_4\) and \(q_6\) parameters and their averaged versions which are widely used in literature. More details can be found here.
import numpy as np
import matplotlib.pyplot as plt
from pyscal3 import System
In this example, we analyse MD configurations, first a set of perfect bcc, fcc and hcp structures and another set with thermal vibrations.
Perfect structures#
We start by creating some reference structures
bcc = System.create.lattice.bcc(lattice_constant=3.147, repetitions=[4,4,4])
fcc = System.create.lattice.fcc(lattice_constant=3.147, repetitions=[4,4,4])
hcp = System.create.lattice.hcp(lattice_constant=3.147, repetitions=[4,4,4])
Next step is calculation of nearest neighbors. There are two ways to calculate neighbors, by using a cutoff distance or by using the voronoi cells. In this example, we will use the cutoff method and provide a cutoff distance for each structure.
Finding the cutoff distance#
The cutoff distance is normally calculated in a such a way that the atoms within the first shell is incorporated in this distance. The :func:pyscal.core.System.calculate_rdf
function can be used to find this cutoff distance.
bccrdf = bcc.calculate.radial_distribution_function()
fccrdf = fcc.calculate.radial_distribution_function()
hcprdf = hcp.calculate.radial_distribution_function()
Now the calculated rdf is plotted
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(11,4))
ax1.plot(bccrdf[1], bccrdf[0])
ax2.plot(fccrdf[1], fccrdf[0])
ax3.plot(hcprdf[1], hcprdf[0])
ax1.set_xlim(0,5)
ax2.set_xlim(0,5)
ax3.set_xlim(0,5)
ax1.set_title('bcc')
ax2.set_title('fcc')
ax3.set_title('hcp')
ax2.set_xlabel("distance")
ax1.axvline(3.6, color='red')
ax2.axvline(2.7, color='red')
ax3.axvline(3.6, color='red')
<matplotlib.lines.Line2D at 0x7fc0fd2e5350>

The selected cutoff distances are marked in red in the above plot. For bcc, since the first two shells are close to each other, for this example, we will take the cutoff in such a way that both shells are included.
Steinhardt’s parameters - cutoff neighbor method#
bcc.find.neighbors(method='cutoff', cutoff=3.6)
fcc.find.neighbors(method='cutoff', cutoff=2.7)
hcp.find.neighbors(method='cutoff', cutoff=3.6)
We have used a cutoff of 3 here, but this is a parameter that has to be tuned. Using a different cutoff for each structure is possible, but it would complicate the method if the system has a mix of structures. Now we can calculate the \(q_4\) and \(q_6\) distributions
bccq = bcc.calculate.steinhardt_parameter([4,6])
fccq = fcc.calculate.steinhardt_parameter([4,6])
hcpq = hcp.calculate.steinhardt_parameter([4,6])
Instead of using the returned values, one can also access them later easily with tab completion
bcc.atoms.steinhardt.generic.q4_norm
array([0.03636965, 0.03636965, 0.03636965, 0.03636965, 0.03636965,
0.03636965, 0.03636965, 0.03636965, 0.03636965, 0.03636965,
0.03636965, 0.03636965, 0.03636965, 0.03636965, 0.03636965,
0.03636965, 0.03636965, 0.03636965, 0.03636965, 0.03636965,
0.03636965, 0.03636965, 0.03636965, 0.03636965, 0.03636965,
0.03636965, 0.03636965, 0.03636965, 0.03636965, 0.03636965,
0.03636965, 0.03636965, 0.03636965, 0.03636965, 0.03636965,
0.03636965, 0.03636965, 0.03636965, 0.03636965, 0.03636965,
0.03636965, 0.03636965, 0.03636965, 0.03636965, 0.03636965,
0.03636965, 0.03636965, 0.03636965, 0.03636965, 0.03636965,
0.03636965, 0.03636965, 0.03636965, 0.03636965, 0.03636965,
0.03636965, 0.03636965, 0.03636965, 0.03636965, 0.03636965,
0.03636965, 0.03636965, 0.03636965, 0.03636965, 0.03636965,
0.03636965, 0.03636965, 0.03636965, 0.03636965, 0.03636965,
0.03636965, 0.03636965, 0.03636965, 0.03636965, 0.03636965,
0.03636965, 0.03636965, 0.03636965, 0.03636965, 0.03636965,
0.03636965, 0.03636965, 0.03636965, 0.03636965, 0.03636965,
0.03636965, 0.03636965, 0.03636965, 0.03636965, 0.03636965,
0.03636965, 0.03636965, 0.03636965, 0.03636965, 0.03636965,
0.03636965, 0.03636965, 0.03636965, 0.03636965, 0.03636965,
0.03636965, 0.03636965, 0.03636965, 0.03636965, 0.03636965,
0.03636965, 0.03636965, 0.03636965, 0.03636965, 0.03636965,
0.03636965, 0.03636965, 0.03636965, 0.03636965, 0.03636965,
0.03636965, 0.03636965, 0.03636965, 0.03636965, 0.03636965,
0.03636965, 0.03636965, 0.03636965, 0.03636965, 0.03636965,
0.03636965, 0.03636965, 0.03636965])
Thats it! Now lets gather the results and plot them.
plt.scatter(bccq[0], bccq[1], s=60, label='bcc', color='#C62828')
plt.scatter(fccq[0], fccq[1], s=60, label='fcc', color='#FFB300')
plt.scatter(hcpq[0], hcpq[1], s=60, label='hcp', color='#388E3C')
plt.xlabel("$q_4$", fontsize=20)
plt.ylabel("$q_6$", fontsize=20)
plt.legend(loc=4, fontsize=15)
<matplotlib.legend.Legend at 0x7fc0fd375790>

Firstly, we can see that Steinhardt parameter values of all the atoms fall on one specific point which is due to the absence of thermal vibrations. Next, all the points are well separated and show good distinction. However, at finite temperatures, the atomic positions are affected by thermal vibrations and hence show a spread in the distribution. We will show the effect of thermal vibrations in the next example.
Structures with thermal vibrations#
Once again, we create the reqd structures Noise can be applied to atomic positions using the noise
keyword as shown below.
bcc = System.create.lattice.bcc(lattice_constant=3.147, repetitions=[4,4,4], noise=0.05)
fcc = System.create.lattice.fcc(lattice_constant=3.147, repetitions=[4,4,4], noise=0.05)
hcp = System.create.lattice.hcp(lattice_constant=3.147, repetitions=[4,4,4], noise=0.05)
cutoff method#
bcc.find.neighbors(method='cutoff', cutoff=3.6)
fcc.find.neighbors(method='cutoff', cutoff=2.7)
hcp.find.neighbors(method='cutoff', cutoff=3.6)
And now, calculate \(q_4\), \(q_6\) parameters
bccq = bcc.calculate.steinhardt_parameter([4,6])
fccq = fcc.calculate.steinhardt_parameter([4,6])
hcpq = hcp.calculate.steinhardt_parameter([4,6])
plt.scatter(fccq[0], fccq[1], s=10, label='fcc', color='#FFB300')
plt.scatter(hcpq[0], hcpq[1], s=10, label='hcp', color='#388E3C')
plt.scatter(bccq[0], bccq[1], s=10, label='bcc', color='#C62828')
plt.xlabel("$q_4$", fontsize=20)
plt.ylabel("$q_6$", fontsize=20)
plt.legend(loc=4, fontsize=15)
<matplotlib.legend.Legend at 0x7fc0fd323950>

We can also visualise the Steinhardt parameter values. For example, the below code plots the \(q_4\) values of the bcc cell.
bcc.show.continuous_property(bcc.atoms.steinhardt.generic.q4_norm)