Calculating coordination numbers

In this example, we will read in a configuration from an MD simulation and then calculate the coordination number distribution.

from pyscal3 import System 
import numpy as np
import matplotlib.pyplot as plt

Read in a file

Lets create a bcc structure

sys = System.create.lattice.bcc(lattice_constant= 4.00, repetitions=[6,6,6])

The above function creates an bcc crystal of 6x6x6 unit cells with a lattice constant of 4.00 along with a simulation box that encloses the particles.

Calculating neighbors

We start by calculating the neighbors of each atom in the system. There are two (main) ways to do this, using a cutoff method and using a voronoi polyhedra method. We will try with both of them. First we try with cutoff system - which has three sub options. We will check each of them in detail.

Cutoff method

Cutoff method takes cutoff distance value and finds all atoms within the cutoff distance of the host atom.

sys.find.neighbors(method='cutoff', cutoff=4.1)

let us try accessing the coordination number of an atom

coordination = [len(sys.atoms.neighbors.index[x]) for x in range(sys.natoms)]

coordination[0]

14

As we would expect for a bcc type lattice, we see that the atom has 14 neighbors (8 in the first shell and 6 in the second). Lets try a more interesting example by reading in a bcc system with thermal vibrations. Thermal vibrations lead to distortion in atomic positions, and hence there will be a distribution of coordination numbers.

sys = System('conf.dump')
sys.find.neighbors(method='cutoff', cutoff=3.6)

We can loop over all atoms and create a histogram of the results

coordination = [len(sys.atoms.neighbors.index[x]) for x in range(sys.natoms)]

Now lets plot and see the results

nos, counts = np.unique(coordination, return_counts=True)
plt.bar(nos, counts, color="#AD1457")
plt.ylabel("density")
plt.xlabel("coordination number")
plt.title("Cutoff method")

Text(0.5, 1.0, 'Cutoff method')

Adaptive cutoff methods

pyscal also has adaptive cutoff methods implemented. These methods remove the restriction on having the same cutoff. A distinct cutoff is selected for each atom during runtime. pyscal uses two distinct algorithms to do this - sann and adaptive. Please check the documentation for a explanation of these algorithms. For the purpose of this example, we will use the adaptive algorithm.

adaptive algorithm

sys.find.neighbors(method='cutoff', cutoff='adaptive', padding=1.5)
coordination = [len(sys.atoms.neighbors.index[x]) for x in range(sys.natoms)]

Now lets plot

nos, counts = np.unique(coordination, return_counts=True)
plt.bar(nos, counts, color="#AD1457")
plt.ylabel("density")
plt.xlabel("coordination number")
plt.title("Cutoff adaptive method")

Text(0.5, 1.0, 'Cutoff adaptive method')

The adaptive method also gives similar results!

Voronoi method

Voronoi method calculates the voronoi polyhedra of all atoms. Any atom that shares a voronoi face area with the host atom are considered neighbors. Voronoi polyhedra is calculated using the Voro++ code. However, you dont need to install this specifically as it is linked to pyscal.

sys.find.neighbors(method='voronoi')

Once again, let us get all atoms and find their coordination

coordination = [len(sys.atoms.neighbors.index[x]) for x in range(sys.natoms)]

And visualise the results

nos, counts = np.unique(coordination, return_counts=True)
plt.bar(nos, counts, color="#AD1457")
plt.ylabel("density")
plt.xlabel("coordination number")
plt.title("Voronoi method")

Text(0.5, 1.0, 'Voronoi method')

Finally..

All methods find the coordination number, and the results are comparable. Cutoff method is very sensitive to the choice of cutoff radius, but voronoi method can slightly overestimate the neighbors due to thermal vibrations.

About small simulation boxes
pyscal repeats the simulation box to create ghost atoms if the input cell is too small. This means that when the neighbors of an atom are extracted, some indices might be repeated.