# Calculating coordination numbers¶

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

import pyscal.core as pc
import numpy as np
import matplotlib.pyplot as plt


The first step is setting up a system. We can create atoms and simulation box using the pyscal.crystal_structures module. Let us start by importing the module.

import pyscal.crystal_structures as pcs

atoms, box = pcs.make_crystal('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. We can then create a System and assign the atoms and box to it.

sys = pc.System()
sys.atoms = atoms
sys.box = box


## Calculating neighbors¶

We start by calculating the neighbors of each atom in the system. There are two 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)


Now lets get all the atoms.

atoms = sys.atoms


let us try accessing the coordination number of an atom

atoms.coordination

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 = pc.System()
sys.find_neighbors(method='cutoff', cutoff=3.6)
atoms = sys.atoms


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

coord = [atom.coordination for atom in atoms]


Now lets plot and see the results

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

Text(0.5, 1.0, 'Cutoff method') 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)
atoms = sys.atoms
coord = [atom.coordination for atom in atoms]


Now lets plot

nos, counts = np.unique(coord, return_counts=True)
plt.ylabel("density")
plt.xlabel("coordination number")

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

atoms = sys.atoms
coord = [atom.coordination for atom in atoms]


And visualise the results

nos, counts = np.unique(coord, return_counts=True)

Text(0.5, 1.0, 'Voronoi method') 