Bond-orientational descriptor

Definition

Bond-order parameters [1] are standard measures of structure in the first coordination shell. Let \(\mathbf{r}_i\) be the position of particle \(i\) and define \(\mathbf{r}_{ij} = \mathbf{r}_j - \mathbf{r}_i\) and \(r_{ij} = |\mathbf{r}_{ij}|\). Then consider the weighted microscopic density around particle \(i\):

\[\rho(\mathbf{r}; i) = \sum_{j=1}^{N_b(i)} w_j \delta(\mathbf{r} - \mathbf{r}_{ij})\]

where \(w_j\) is a particle-dependent weight and the sum involves a set of \(N_b(i)\) particles, which defines the coordination shell of interest for particle \(i\).

We project the microscopic density on a unit-radius sphere, that is, \(\hat{\rho}(\hat{\mathbf{r}}; i) = \sum_{j=1}^{N_b(i)} w_j \delta(\mathbf{r} - \hat{\mathbf{r}}_{ij})\), where \(\hat{\mathbf{r}} = \mathbf{r} / |\mathbf{r}|\) and similarly \(\hat{\mathbf{r}}_{ij} = \mathbf{r}_{ij}/|\mathbf{r}_{ij}|\). Expanding in spherical harmonics yields

\[\hat{\rho}(\hat{\mathbf{r}}; i) = \sum_{l=0}^\infty \sum_{m=-l}^l c_{l m}(i) Y_{l m}(\hat{\mathbf{r}}) ,\]

with coefficients

\[c_{l m}(i) = \int d\mathbf{r} \rho(\mathbf{r}; i) Y_{l m}(\hat{\mathbf{r}}) .\]

In the conventional bond-order analysis, one sets the weights \(w_j\) to unity and considers the normalized complex coefficients,

\[\begin{split}\begin{align} q_{lm}(i) & = \frac{1}{N_b(i)} \int d\mathbf{r} \rho(\mathbf{r}; i) Y_{l m}(\hat{\mathbf{r}}) \nonumber \\ & = \frac{1}{N_b(i)} \sum_{j=1}^{N_b(i)} Y_{l m}(\hat{\mathbf{r}}_{ij}) . \end{align}\end{split}\]

The rotational invariants,

\[Q_{l}(i) = \sqrt{ \frac{4\pi}{2l + 1}\sum_{m=-l}^l |q_{lm}(i)|^2 },\]

provide a detailed structural description of the local environment around particle \(i\).

We then consider \(Q_l(i)\) for a sequence of orders \(\{ l_n \} = \{ l_\mathrm{min}, \dots, l_\mathrm{max} \}\). The resulting feature vector for particle \(i\) is given by

\[X^\mathrm{BO}(i) = (\: Q_{l_\mathrm{min}}(i) \;\; \dots \;\; Q_{l_\mathrm{max}}(i) \:) .\]

Setup

Instantiating this descriptor on a Trajectory can be done as follows:

from partycls import Trajectory
from partycls.descriptors import BondOrientationalDescriptor

traj = Trajectory("trajectory.xyz")
D = BondOrientationalDescriptor(traj)

The constructor takes the following parameters:

BondOrientationalDescriptor.__init__(trajectory, lmin=1, lmax=8, orders=None, accept_nans=True, verbose=False)[source]
Parameters

  • trajectory (Trajectory) – Trajectory on which the structural descriptor will be computed.

  • lmin (int, default: 1) – Minimum order \(l_\mathrm{min}\). This sets the lower bound of the grid \(\{ l_n \}\).

  • lmax (int, default: 8) – Maximum order \(l_\mathrm{max}\). This sets the upper bound of the grid \(\{ l_n \}\). For numerical reasons, \(l_\mathrm{max}\) cannot be larger than 16.

  • orders (list, default: None) – Sequence \(\{l_n\}\) of specific orders to compute, e.g. orders=[4,6]. This has the priority over lmin and lmax.

  • accept_nans (bool, default: True) – If False, discard any row from the array of features that contains a NaN element. If True, keep NaN elements in the array of features.

  • verbose (bool, default: False) – Show progress information and warnings about the computation of the descriptor when verbose is True, and remain silent when verbose is False.

Hint

The alias SteinhardtDescriptor can be used in place of BondOrientationalDescriptor.

Requirements

The computation of this descriptor relies on:

  • Lists of nearest neighbors. These can either be read from the input trajectory file, computed in the Trajectory, or computed from inside the descriptor using a default method.

Demonstration

We consider an input trajectory file trajectory.xyz in XYZ format that contains two particle types "A" and "B". We compute the lists of nearest neighbors using the fixed-cutoffs method:

from partycls import Trajectory

# open the trajectory
traj = Trajectory("trajectory.xyz")

# compute the neighbors using pre-computed cuttofs
traj.nearest_neighbors_cuttofs = [1.45, 1.35, 1.35, 1.25]
traj.compute_nearest_neighbors(method='fixed')
nearest_neighbors = traj.get_property("nearest_neighbors")

# print the first three neighbors lists for the first trajectory frame
print("neighbors:\n",nearest_neighbors[0][0:3])
Output:
neighbors:
 [list([16, 113, 171, 241, 258, 276, 322, 323, 332, 425, 767, 801, 901, 980])
  list([14, 241, 337, 447, 448, 481, 496, 502, 536, 574, 706, 860, 951])
  list([123, 230, 270, 354, 500, 578, 608, 636, 639, 640, 796, 799, 810, 826, 874, 913])]

We now instantiate a BondOrientationalDescriptor on this trajectory and restrict the analysis to type-B particles only. We set set the grid of orders \(\{l_n\} = \{2,4,6,8\}\):

from partycls.descriptors import BondOrientationalDescriptor

# instantiation
D = BondOrientationalDescriptor(traj, orders=[2,4,6,8])

# print the grid of orders
print("grid:\n", D.grid)

# restrict the analysis to type-B particles
D.add_filter("species == 'B'", group=0)

# compute the descriptor's data matrix
X = D.compute()

# print the first three feature vectors
print("feature vectors:\n", X[0:3])
Output:
grid:
 [2 4 6 8]
feature vectors:
 [[0.06498973 0.10586717 0.46374576 0.22207796]
  [0.12762569 0.09640384 0.49318559 0.29457554]
  [0.08327171 0.11151433 0.37917788 0.17902556]]

  • grid shows the grid of orders \(\{ l_n \}\).

  • feature vectors shows the first three feature vectors \(X^\mathrm{BO}(1)\), \(X^\mathrm{BO}(2)\) and \(X^\mathrm{BO}(3)\) corresponding to the grid.

References

1

Paul J. Steinhardt, David R. Nelson, and Marco Ronchetti. Bond-orientational order in liquids and glasses. Phys. Rev. B, 28(2):784–805, 1983. doi:10.1103/PhysRevB.28.784.