Smoothed bond-orientational descriptor

Important

See Bond-orientational descriptor first.

Definition

This is a smooth version of the Bond-orientational descriptor, in which the coefficients \(q_{lm}(i)\) are multiplied by a weighting function \(f(r)\) that depends on the radial distance \(r\) between the central particle \(i\) and other surrounding particles \(j\), where \(j\) can be any particle in the system (i.e. not necessarily a nearest neighbors of \(i\)).

The smoothed complex coefficients are given by

\[q_{lm}^{S}(i) = \frac{1}{Z(i)} \sum_{j=1}^{N} f({r}_{ij}) Y_{lm}(\hat{\mathbf{r}}_{ij}) ,\]

where \(Z(i)=\sum_{j=1}^{N} f({r}_{ij})\) is a normalization constant and the superscript \(S\) indicates the smooth nature of the descriptor. We use

\[f(r_{ij}) = \exp \left[- (r_{ij} / r_{\alpha\beta}^c)^\gamma \right] H(R_{\alpha\beta}^c - r_{ij}) ,\]

where \(r_{\alpha\beta}^c\) is the first minimum of the corresponding partial radial distribution function for the pair \((i,j)\) and \(\gamma\) is an integer. Also, \(H\) is the Heaviside step function, which ensures, for efficiency reasons, that the descriptor only has contributions from particles within a distance \(R_{\alpha\beta}^c = \xi \times r_{\alpha\beta}^c\) from the central one, where \(\xi > 1\) is a scaling factor.

The rotational invariants are defined similarly to the Bond-orientational descriptor.

We then consider \(Q_l^S(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{SBO}(i) = (\: Q_{l_\mathrm{min}}^S(i) \;\; \dots \;\; Q_{l_\mathrm{max}}^S(i) \:) .\]

Setup

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

from partycls import Trajectory
from partycls.descriptors import SmoothedBondOrientationalDescriptor

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

The constructor takes the following parameters:

SmoothedBondOrientationalDescriptor.__init__(trajectory, lmin=1, lmax=8, orders=None, cutoff_enlargement=1.3, exponent=8, 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.
cutoff_enlargement (float, default: 1.3) – Scaling factor \(\xi\) for the nearest neighbors cutoffs \(r_{\alpha\beta}^c\) to consider neighbors \(j\) a distance \(R_{\alpha\beta}^c = \xi \times r_{\alpha\beta}^c\) away from the central particle \(i\).
exponent (int, default: 8) – Exponent \(\gamma\) in the smoothing function \(f(r_{ij})\).
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.

Requirements

The computation of this descriptor relies on:

Nearest neighbors cutoffs. These can either be set in the Trajectory prior to the computation of the descriptor, 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 can either compute or set the nearest neighbors cutoffs \(\{ r_{\alpha\beta}^c \}\) for the smoothing directly in Trajectory:

from partycls import Trajectory

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

# compute the nearest neighbors cutoffs
traj.compute_nearest_neighbors_cutoffs(dr=0.1)
print("computed cuttofs\n", traj.nearest_neighbors_cutoffs)

# set the nearest neighbors cuttofs
traj.nearest_neighbors_cutoffs = [1.45, 1.35, 1.35, 1.25]
print("manually set cuttofs\n", traj.nearest_neighbors_cutoffs)

Output:

computed cutoffs:
 [1.4500000000000004, 1.3500000000000003, 1.3500000000000003, 1.2500000000000002]
manually set cutoffs:
 [1.45, 1.35, 1.35, 1.25]

Note

If not computed in Trajectory or manually set, the cutoffs \(\{ r_{\alpha\beta}^c \}\) will be computed from inside the descriptor.

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

from partycls.descriptors import SmoothedBondOrientationalDescriptor

# instantiation
D = SmoothedBondOrientationalDescriptor(traj,
                                        orders=[2,4,6,8],
                                        cutoff_enlargement=1.3,
                                        exponent=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.03156284 0.11095233 0.4112718  0.23829355]
  [0.0698711  0.08107918 0.47678647 0.26671868]
  [0.06221017 0.09806095 0.39152213 0.16630718]]

grid shows the grid of orders \(\{ l_n \}\).
feature vectors shows the first three feature vectors \(X^\mathrm{SBO}(1)\), \(X^\mathrm{SBO}(2)\) and \(X^\mathrm{SBO}(3)\) corresponding to the grid.