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\):
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
with coefficients
In the conventional bond-order analysis, one sets the weights \(w_j\) to unity and considers the normalized complex coefficients,
The rotational invariants,
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
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 overlmin
andlmax
.accept_nans (bool, default: True) – If
False
, discard any row from the array of features that contains a NaN element. IfTrue
, 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 isFalse
.
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])
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])
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.