An atomistic fingerprint algorithm for learning ab initio molecular force fields

Molecular Dynamics (MD) simulations have been widely used for studying atomistic systems, e.g., proteins and catalysts, due to their ability to precisely capture transient events and to predict macroscopic properties from microscopic details.^1,2 In its most prevalent implementation, the trajectory of an atomistic system is integrated in time according to the Newton’s law of motion using forces calculated as the negative gradient of a Hamiltonian, whose functional form and parameters are collectively referred to as a force field.^3–5 Traditionally, the pairwise and many-body terms that comprise a force field are derived empirically by fitting to quantum mechanical (QM) calculations and experimental data.

Three properties directly relate to the applicability of a force field: accuracy, transferability, and complexity.^6,7 Over the years, a large number of force fields have been developed, each carrying a particular emphasis over these three properties. However, the combinatorial complexity of atomistic systems can easily outpace force field development efforts, the difficulty of which explodes following the curse of dimensionality.⁸ A deceptively simple system that can demonstrate the situation is water, a triatomic molecule with a well-characterized molecular structure. In fact, all common water models, such as SPC-E, TIP3P, and TIP4P, have only succeeded in reproducing a small number of structural and dynamical properties of water due to the difficulty in modeling strong intermolecular many-body effects such as hydrogen bonding and polarization.^9,10

In lieu of a force field, quantum mechanical (QM) calculations can be employed straightforwardly to drive molecular dynamics simulations. The method achieves significantly better accuracy and transferability by solving for the electronic structure of the system. However, the computational complexity of QM methods is at least cubic in the number of electrons, and consequently the time and length scales accessible by QM-driven molecular dynamics are severely constrained.

Assuming that there is smoothness in the potential energy surface of the atomistic system, one possible strategy to accelerate QM-driven molecular dynamics is to use QM calculations on only a subset of the time steps and to interpolate for similar atomic configurations.^11–14 A schematic overview of the process is given in Fig. 1, which is enabled by the recent development of high-dimensional nonlinear statistical learning and regression techniques such as Gaussian process regression (GPR)¹⁵ and artificial neural networks.¹⁶ 

This paper focuses on a particular aspect of the machine-learning-driven molecular computation pipeline, i.e., fingerprint algorithms, whose importance arises naturally from the aforementioned regression protocol. A fingerprint is an encoding of an atomistic configuration that can facilitate regression tasks such as similarity comparison across structures consisting of variable numbers of atoms and elements. As has been pointed out previously,¹² a good fingerprint should possess the following properties:

1. It can be encoded as a fixed-length vector so as to facilitate regression (particularly for artificial neural networks).

2. It is complete, i.e., different atomistic neighborhood configurations lead to different fingerprints and vice versa, and the “distance” between the fingerprints should be proportional to the intrinsic difference between the atomistic neighborhood configurations.

3. It is continuous with regard to atomistic coordinates, and the change in fingerprint should be approximately proportional to the structural variation as characterized by, for example, some internal coordinates.

4. It is invariant under permutation, rotation, and translation.

5. It is computationally feasible and straightforward to implement.

Before we proceed to the details of our work, we will first briefly review several fingerprints that are closely related to our work, i.e., the Smooth Overlap of Atomic Positions (SOAP) kernel,¹² the Coulomb matrix,¹⁷ and the Graph Approximated Energy (GRAPE) kernel.¹⁸ 

Smooth Overlap of Atomic Positions (SOAP): The SOAP kernel is built on the idea of representing atomistic neighborhoods as smoothed density fields using Gaussian kernels each centered at a neighbor atom. Similarity is measured as the inner product between density fields, while rotational invariance is achieved by integrating over all possible 3D rotations, which can be performed analytically using the power spectrum of the density field. In fact, our fingerprint algorithm is inspired by this idea of treating atoms as smoothed density fields. However, we take a different approach to endorse the fingerprint with rotational invariance and use the Euclidean distance instead the of inner product as a distance metric.

Coulomb matrix: The practice of using graphs to represent atomistic neighbor configurations was first implied by the Coulomb matrix, and later further formulated in the GRAPE kernel, where the diffusion distance was proposed as a similarity measure between different local chemical environments.¹⁹ The idea is to construct an undirected, unlabeled graph G = (V, E) with atoms serving as the vertices and pairwise interactions weighting the edges. For example, the Coulomb matrix can be treated as a physically inspired Laplacian matrix,²⁰ 

where the degree matrix D encodes a polynomial fit of atomic energies to the nuclear charge, while the adjacency matrix A corresponds to the Coulombic interactions between all pairs of atoms. Due to the use of only relative positions between atoms in the adjacency matrix, the Coulomb matrix is automatically invariant under translation and rotation. However, the matrix itself is not invariant under permutation, as swapping the order of two atoms will result in an exchange of the corresponding columns and the rows. To address this, the sorted list of eigenvalues of the Coulomb matrix can be used instead as a feature vector, while an ℓ_p norm can be used as a distance metric. In practice, due to the fact that the number of neighbor atoms may change, the shorter eigenvalue list is padded with zeros in a distance computation.

Graph Approximated Energy (GRAPE): The GRAPE kernel evaluates the simultaneous random walks on the direct product of the two graphs representing two atomistic neighborhood configurations. Permutational invariance is achieved by choosing a uniform starting and stopping distribution across nodes of both graphs. However, the cost of distance computation between two graphs scales as $O(N2)$ with a one-time per-graph diagonalization cost of $O(N3)$⁠.

In Secs. II–V, we present our new fingerprint algorithm, namely, the Density-Encoded Canonically Aligned Fingerprint (DECAF). The paper is organized as follows: in Sec. II, we introduce a robust algorithm that can determine canonical coordinate frames for obtaining symmetry-invariant projections; in Sec. III, we present numerical recipes to use smoothed atomistic density fields as a fingerprint for molecular configuration; in Sec. IV, we demonstrate the capability of the fingerprint via examples involving the regression of atomistic potential energy surfaces; in Sec. V, we discuss the connection between our algorithm and previously proposed ones; and we conclude with a discussion in Sec. VI.

To improve model generalization while minimizing data redundancy, a fingerprint algorithm should be able to recognize atomistic structures that differ only by a rigid-body transformation or a permutation of atoms of the same element and to extract feature vectors invariant under these transformations. As summarized in Table I, a variety of strategies have been successfully employed by common fingerprint algorithms to achieve rotational invariance.

However, these approaches do not provide a means for the acquisition of vector-valued quantities in a rotational invariant form. One approach is to only acquire and interpolate the potential energy, a scalar quantity, and then take the derivative of the regression model. This approach, however, triggers the need for methods to decompose the total energy among the atoms, which is a property of the entire system rather than individual atoms.²¹ 

Another approach proposed by Li et al.^13,22 is to project vector quantities onto a potentially overcomplete set of non-orthogonal basis vectors obtained from a weighted sum of the atomic coordinate vectors,

However, the approach may suffer from robustness issues. For example, all of the V_k generated with different p_k will point in the same direction if the radial distances of the atoms are all equal. Further, the configuration with 4 atoms at (r cos ε, r sin ε), (0, r), (−r, 0), and (0, −r) leads to

Thus, if ε gets close to zero, V_k will always point toward either (0, 1) or (0, −1), even if the vector quantity of interest may point in other directions.

Here, we present a robust kernel principal components analysis (PCA)-inspired algorithm for the explicit determination of a canonical coordinate frame, within which the projection of the atomistic neighborhood is invariant under rigid-body rotation. Furthermore, the canonical coordinate frame can be directly used to capture vector-valued quantities in a rotational-invariant form. Given N atoms with position $x1,…,xN∈Rd$⁠, we first formulate the L_p PCA algorithm as an optimization problem where we seek a unit vector w* that maximizes the sum of the projections,

where r_i = ∥x_i∥ is the distance from the origin to atom i and e_i = x_i/r_i is the unit vector pointing toward atom i, respectively. The optimization process can only uniquely determine the orientation of a projection vector up to a line because $wTe≡−wTe$⁠. As a consequence, further heuristics are needed to identify a specific direction for the PCA vectors.

To overcome this difficulty, we generalize the $rip$ term into a weight function g(r) of radial distance and the $wTeip$ term into a bivariate kernel function κ(w, e) between two vectors. We then attempt to seek a unit vector w* that minimizes the kernel summation,

We will hereafter refer to this formulation as kernel minisum.

In particular, we have found a square angle (SA) kernel and an exponentiated cosine (EC) kernel that perform well in practice,

As shown in Fig. 2, both kernels are minimal when w and e are parallel and monotonically reach maximum when w and e are antiparallel. Intuitively, optimizing the minisum objective function generated by the SA kernel will yield a vector that, loosely speaking, bisects the sector occupied by the atoms. The EC kernel exhibits very similar behavior but leads to a smoother objective function. As shown in Fig. 2, this allows for the determination of a projection vector without ambiguity, even if the atom configuration contains perfect symmetry.

A major advantage of the kernel minisum approach versus L^p norm-based PCA lies in its (1) robustness in the presence of structural symmetry and (2) continuity of the resulting principal axes with respect to angular movement of the input data. As shown in Fig. 3, kernel minisum is particularly suitable for atomistic systems where strong symmetries are common and the continuity against angular movement is desired. The minisum framework can also be used with other customized kernels to suit for the characteristics of specific application scenarios.

The optimization problem can be solved very efficiently using a gradient descent algorithm as detailed below.

Square Angle: The objective function of the minisum problem using the square angle (SA) kernel is

The gradient of A_SA with respect to w is

Note that $arccos(wTei)1−(wTei)2$ is singular when w∥e_i. This can be treated numerically by replacing the removable singularities at w^Te_i = 1 with the left-limit $limwTei→1−arccos(wTei)1−(wTei)2=1$⁠, while truncating the gradient at a finite threshold near the poles at w^Te_i = −1.

A local minimum can be iteratively searched for with gradient descent while renormalizing w after each iteration. Moreover, due to the locally quadratic nature of the objective function, we have found that the Barzilai-Borwein algorithm²³ can significantly accelerate the convergence at a minimal cost. The algorithm is presented in Algorithm 1.

Exponentiated Cosine: The objective function of the minisum problem using the exponentiated cosine (EC) kernel is

The gradient of A_EC with respect to w is

The gradient contains no singularity. However, it is not always locally quadratic or convex. This can cause the Barzilai-Borwein algorithm to generate negative step sizes and consequently divert the search towards a maximum. Luckily, this can be easily overcome by always using the absolute value of the step size generated by the Barzilai-Borwein algorithm. Such enforcement prevents the minimization algorithm from going uphill. The complete algorithm is given in Algorithm 2.

As shown in Table II, both Algorithms 1 and 2 converge quickly and consistently across a wide range of representative point configurations commonly found in molecular systems. However, the gradient descent method can only find local optima. Thus, multiple trials should be performed using different initial guesses to ensure that a global minimum can be located.

In 3D, a complete set of orthogonal bases can be found greedily using the protocol as described in Algorithm 3. Specifically, we use the globally optimal solution of the minisum optimization problem as the first basis b_α, and the constrained optimal solution in a plane orthogonal to b_α as the second basis b_β. Special care must be taken for determining the third basis b_γ, as the only degree of freedom now is its sign due to the orthogonality constraint. The straightforward approach of choosing the direction that gives the smaller objective function value may fail, for example, when the system contains improper rotational symmetry. In that case, b_α and b_β are interchangeable and both are perpendicular to the rotation-reflection axis. As a result, the two candidates of b_γ will both align with the rotation-reflection axis and are thus indistinguishable by kernel minisum. However, the projection of the atoms onto the two seemingly equivalent coordinate frames is not identical but rather mirroring images of each other. Fortunately, this can be addressed by choosing the direction of the half-space, as created by the plane b_α-b_α, that yields the smaller kernel objective function between the bisector of b_α and b_β versus the points which lie in that half-space. This rule can also handle general situations with/without symmetry.

It is difficult to prove global uniqueness of the kernel minisum solution given the non-convex nature of the exponentiated cosine and square angle kernels. In fact, it seems that the only kernel that allows analytical proof of solution uniqueness is κ_COM(w, e) ≐ −w^Te, whose solution simply corresponds to the weighted center of mass of the neighbor atoms. However, this simple kernel is not robust against reflectional and rotational symmetry. Luckily, the rare cases where two global optimal solutions do coexist can be safely captured by the repeated searching procedure starting from different seeds. Thus, a fingerprint can be extracted using each of the resulting coordinate frame. This may mildly increase the size of the training set, which could even be advantageous when training data are scarce.

The local density field ρ_s(r) around a point s, as shown in Fig. 4, is formulated as a superimposition of smoothing kernel functions each centered at a neighbor atom i = 1, 2, …, N with relative displacement x_i with regard to s and within a cutoff distance R_c,

This density field, as has been pointed out previously,¹² may be used as a fingerprint of the local atomistic environment. Here, we assume that the smoothing kernel W_d(r) takes the form of a stationary Gaussian $σ−1exp−12∥r∥2/σ2$⁠. We also assume that the density scaling function W_c(r), which ensures the continuity of the density field when atoms enter or exit the cutoff distance, is a bell-shaped function with compact support. Further discussion on both W_c(r) and W_d(r) can be found in Secs. III B and III C, respectively.

To achieve rotational invariance, we project the atom coordinates into the canonical coordinate frame R ≐ [b_α, b_β, b_γ] as determined by the kernel minisum algorithm, when generating the density field,

Depending on the specific application, s may not necessarily overlap with any of the x_i. Scalar properties can be acquired directly from the target atom, while vector-valued properties, such as force, can be acquired and interpolated in the local orthogonal coordinates as $ỹ$ = R^Ty.

We define the distance between two density fields ρ_i and ρ_j as a weighted volume integral of their pointwise difference,

The weight function $w$(r) provides additional flexibility for emphasizing particular regions of the atomistic neighborhood. It could play an important role for fitting properties with steep gradients, e.g., the repulsive part of the Lennard-Jones potential.

We now introduce an optimal quadrature rule to approximate the integral in Eq. (19) in a computationally tractable manner. A quadrature rule is a numerical recipe in the form $∫f(x)=∑i=0Nwif(xi)$⁠, which numerically approximates a definite integral using only discrete evaluations of the integrand. To determine the quadrature nodes and weights, we decompose the volume integral in Eq. (19) into a surface integral over spherical shells and a 1D integral along the radial direction,

The surface integral can be optimally approximated using the Lebedev quadrature rule,²⁴ 

where β_m, q_m ≐ (x_m, y_m, z_m), and N_a are the weights, positional unit vectors, and number of the Lebedev nodes, respectively. The radial integral fits well into the generalized Laguerre-Gauss quadrature formula with weight function r²e^−r,²⁵ 

where α_n, r_n, and N_r are the weights, coordinates, and number of the Laguerre nodes, respectively. Combining Eqs. (20)–(22), a composite quadrature rule can be generated consisting of several spherical layers of nodes. As shown in Fig. 5, the radial coordinates of the quadrature nodes are determined by the Laguerre quadrature nodes, while the angular positions are determined by the Lebedev quadrature nodes, respectively. This composite quadrature formula translates the 3D volume integral into a summation over discrete grid points,

Using the right hand side of Eq. (23) to replace the integral in Eq. (19) and using the multi-index notation k = (n, m), 1 ≤ n ≤ N_r, 1 ≤ m ≤ N_a(n) to enumerate over the quadrature nodes located at r_k = r_n ⋅q_m with weights $wk=4παnβmw(rn⋅qm)ern$⁠, we obtain the final discretized formula for computing the distance between the fingerprints,

For quick reference, we tabulated in the  Appendix the values for r_n and α_n in the Laguerre quadrature of up to 6 points and the values for q_m, β_m in the Lebedev quadrature of up to 50 points.

In addition, the quadrature nodes could be radially scaled such that the outer most nodes lie at a radius R* within some cutoff distance R_c. This allows us to fit a Laguerre quadrature of any order within an arbitrary cutoff distance. The scaled quadrature rule is given by

Since the scaling is simply constant among all nodes, it can be safely omitted in many regression tasks where only the relative proximity between the fingerprints are of significance.

In this section, we examine two radial weight functions that can be used to fine-tune the density field fingerprint: the density scaling function W_c(r) as appears in Eq. (17) and the weight of integral w(r) as appears in Eq. (19).

Driven by the interest of reducing computational cost, we would like to use a cutoff distance to select atoms involved in constructing the density field. However, it is important to ensure that atoms will enter and exit the neighborhood smoothly. This naturally requests that the contribution of an atom to be zero outside of the cutoff and to increase continuously and smoothly when the atom approaches entrance. Correspondingly, W_c(r) should (1) become unity at the origin, (2) smoothly approach zero at the cutoff distance, and (3) be twice-differentiable to ensure the differentiability of regression models based on the fingerprint. Candidates satisfying the above conditions include, for example, a tent-like kernel

and a bell-shaped polynomial kernel with compact support

as detailed in the  Appendix.

The approximation of the radial integral using a Laguerre quadrature requires that the integrand, i.e., the pointwise difference between the atomistic density fields, decays sufficiently fast beyond the outermost quadrature nodes in order to achieve acceptable convergence. In addition, the steeply repulsive yet flatly attractive interatomic short-range interactions prompt that the sensitivity of fingerprint be adjusted correspondingly in order to avoid numerical difficulties in training a regression model. The weight of the integral, $w$(r), provides a convenient means for achieving the purpose. Different from W_c(r), $w$(r) should instead satisfy the following conditions: (1) should be normalized such that ∫$w$(r)dV (r) = 1; (2) should decay sufficiently fast, but not necessarily to 0, beyond the outer most quadrature node; and (3) should be sufficiently smooth beyond the outermost quadrature node. Candidates for $w$(r) include the tent-like kernel and the bell-shaped kernel for W_c(r), albeit with a different normalization factor. The Laplacian kernel

with a properly sized length scale l also appears to be a candidate due to its similarity with e^−r part of the weight function of the Laguerre quadrature. Note that the constant kernel

may also be a choice as long as the density field already decays fast enough due to the density scaling function W_c(r).

In Fig. 6, we demonstrate the effect of the density scaling function and the weight of integral on the distance matrices between fingerprints obtained from a pair of atoms. The comparison between panels (a) and (b) shows that a bell-shaped integration weight allows the distance between fingerprints to change more rapidly when the atoms are closer but more slowly when the atoms are farther apart. The visible discontinuity in the second row clearly demonstrates the importance of the damping function when only atoms within a finite cutoff distance are used to compute the fingerprint.

We further examine the impact of the weight functions on the performance of Gaussian process regression (GPR) using the fingerprint of the interatomic force of a minimal system containing two nitrogen atoms. Despite the simplicity of the system, this case is of fundamental importance because of its ubiquity and because the fast-growing repulsive regime of the Lennard-Jones potential could cause difficulty as a small change in the system configuration can trigger large changes in the regression target function. In Fig. 7, we compare the performance among the combination of four weights of the integral and two density scaling functions. The initial training set consists of two samples collected at r_N–N = 1.0 and 6.0. The regression is then refined using a greedy strategy that consecutively learns the point with the largest posterior variance until the largest uncertainty, defined as twice the posterior standard deviation, is less than 0.1 eV/Å. The active learning scheme is able to deliver a GPR model, for each and every combination of the weight functions, that closely fits the target function. However, the numbers of refinement iterations and the achieved accuracy do vary. Therefore, it is important to evaluate and choose the weight functions in the context of specific application scenarios.

Despite the formal convergence of the composite quadrature in DECAF, a cost of $O(NM)$ distance calculations and kernel evaluations is needed to sample a density field generated by N atoms using M quadrature nodes. A less prominent cost is associated with the L² distance calculation, which comes at a cost of $O(M)$ floating point operations. Thus, in practice, it is often desirable to use as few nodes as possible to capture only information of the density field within a certain band limit.²⁶ Accordingly, the integral cutoff R_c, the number of quadrature nodes, and the width of the density kernel need to be tuned to obtain an optimal balance between resolution and computational speed.

When designing the composite quadrature rule, we chose the Laguerre quadrature for the radial direction because its nodes are denser near the origin but sparser farther away. This is consistent with our physical intuition that the near field generally has a stronger influence than the far field in an atomistic neighborhood. For example, the van der Waals potential grows rapidly when atoms are in direct contact but flattens out of the first coordinate shell. Accordingly, it may be possible for us to use sparser outer-layer grids to reduce the total number of quadrature nodes, while still keeping enough nodes in the inner layers to maintain the sensitivity of the quadrature toward close neighbors. Cooperatively, we can also use non-stationary Gaussian density kernels whose width is dependent on the distance from the atom to the origin. In this way, the sparser nodes should still sufficiently sample the smoother far field. Wider kernels at remote atoms also reduce the total difference between the far fields of two fingerprints in a statistical sense. Thus, the contribution of the far field in the integral can be effectively tuned even though the weights on the quadrature nodes remain the same.

In Fig. 8, we demonstrate how a variable-resolution quadrature can be combined with a widening smoothing density kernel to simultaneously reduce the computational complexity and preserve the quality of the fingerprint. In column A, a dense grid is used to sample density fields generated by a wide smoothing length. By examining the distance matrices of fingerprints sampled during bond stretching and angular stretching movements, we note that the radial similarity decreases monotonically while the angular similarity changes nearly constantly. In column B, the number of quadrature nodes is kept the same, but the smoothing length is reduced as an attempt to increase the fingerprint sensitivity. Better response in the near field of the radial direction is obtained, but the linearity in the far field in the angular direction is compromised. In column C, the fingerprint performs even worse due to the combination of a sparser quadrature grid and a small smoothing length. In column D, the performance is recovered because we let the smoothing length parameter σ of the Gaussian density kernels W_d(r) depend on the distance from the origin to each atom and simultaneously adjust the quadrature node density according to this pattern.

Regression tasks throughout this work are performed using Gaussian process regression (GPR), a nonlinear kernel method, that treats training data points as discrete observations from an instantiation of a Gaussian process. Predictions are made using the posterior mean and variance of the joint distribution between the test data and the training data. One particular interesting property about the Gaussian process is that the posterior variance may be interpreted as a measure of prediction uncertainty, which can be exploited to design active learning algorithms for sampling expensive functions. The actual computation used our own software implementation which was made publicly available on Zenodo.²⁷ We use the square exponential covariance kernel to compute the covariance, i.e., correlation, between the samples,

where x and x′ are the DECAF fingerprints, and d(x, x′) is the distance between norms as computed by Eq. (24) or Eq. (26). The kernel is stationary, meaning that the covariance depends only on the relative distance between two samples but not their absolute position. The training process searches for the hyperparameters, i.e., the output variance σ and the length scale l, that maximizes the likelihood of the training data. A detailed tutorial on GPR can be found in Ref. 15. An illustration on the complete workflow of using the density field fingerprint to perform regression tasks is given in Fig. 9.

First, we attempt to fit the potential energy surface of a protonated water dimer system, in a head-to-head configuration, as a function of the oxygen-oxygen distance r_O–O and the dihedral angle φ between the two planes each formed by a water molecule. As shown in Fig. 10(a), the system contains an improperly rotational symmetry, which we wish to capture with the kernel minisum algorithm. A GPR model was seeded with 8 training points corresponding to the combinations of r_O–O = 2.2, 2.4, 2.6, 2.8 and φ = 0, π/2. Subsequently, an active learning protocol was used to greedily absorb points with the highest uncertainty into the training set. Despite that we restricted all training data to be within the subdomain φ ≤ π/2, as shown by Figs. 10(b) and 10(c), we are able to accurately reproduce the target function over the entire parameter space after a few active learning steps.

The DECAF fingerprint used here is constructed with 3 spherical layers within a cutoff distance R_c of 6.0 Å, each consisting of 14, 26, and 38 Lebedev quadrature nodes, respectively. The weight of integral was chosen as $w$(r) = W^6,4(1 − r/R_c), where W^6,4 is the bell-shaped polynomial as defined in the  Appendix, Eq. (A1). The density scaling function $Wc(r)=(1−∥r∥/Rc)3$⁠, where r is the vector from the atom to the fingerprint center, is the tent-like kernel as defined in Eq. (27) with t = 3. The density kernel that sits on the oxygen atoms assumes the form of a non-stationary Gaussian as discussed in Sec. III C: $WdO(r,r′)=σO(r)−1exp−12∥r′∥2/σO(r)2,σO(r)=1.5+0.25∥r∥$ with r being the vector from the atom to the fingerprint center and r′ being the vector from the atom to the quadrature node. The density kernel for the hydrogen atoms has a different weight and width to ensure discriminability: $WdH(r,r′)=0.75σH(r)−1exp−12∥r′∥2/σH(r)2,σH(r)=0.9+0.15∥r∥$⁠.

Next, we demonstrate the usability of fingerprint for fitting vector-valued quantities by performing geometry optimization and vibrational analysis on a single water molecule. The process involves the simultaneous regression of: (1) energy, a molecular scalar quantity; (2) force, a per-atom vector quantity; and (3) dipole, a molecular vector quantity. Correspondingly, we performed GPR of energy and dipole using fingerprints extracted from the center of mass of the molecule and GPR of force using fingerprints extracted from each atom. Each component of the vector properties is modeled independently as a scalar Gaussian process. The training set consists of 45 configurations uniformly covering the range r_O–H = 0.93, 0.95, 0.97, 0.99, 1.05 Å and θ_H–O–H = 101°, 105.5°, 111°. As shown in Table III, the GPR model can successfully drive calculations of the infrared spectrum of the molecule from randomly perturbed initial structures in arbitrary orientation. The fingerprint configuration is the same as in Sec. IV B.

As shown in Fig. 11, here we attempt to fit for the forces felt by the atoms in a benzene molecule along the MD Trajectories as obtained from a sibling database of QM7.^28,29 The density kernel for the carbon atoms assumes the same functional form with that of the oxygen atoms, but uses a different smoothing length function σ_C(r) = 1.2 + 0.2∥r∥. The rest of the parameters are inherited from the previous examples. The training configurations were chosen adaptively in an iterative process using the sum of the GPR posterior variance and the prediction error as the acquisition function.

In Fig. 12, we compare the ability to distinguish atomistic configurations of our fingerprint with SOAP and the Coulomb matrix. Our work is inspired by the SOAP descriptor,¹² which proposes the use of smoothed densities to represent atomistic neighborhoods. However, instead of converting the density field into the frequency domain using spherical harmonics, we perform density field sampling and comparison directly in the real space. This is enabled due to the availability of the canonical coordinate frame as computed through the kernel minisum optimization. We have mainly used the L₂ norm to compute the distance between atomistic neighborhoods. However, our fingerprint exhibits very similar behavior to SOAP when used together with an inner product formula

as demonstrated in Fig. 12(a). Thus, our fingerprint could be used in conjunction with a wide variety of covariance functions based on either the Euclidean distance or the inner product similarity.

At first sight, DECAF is very different from the Coulomb matrix fingerprint and GRAPE, which are both graph-based algorithms.^17,18 However, instead of trying to capture the overall density field, if we measure the contribution from each individual atom on the quadrature nodes at z₁, z₂, …, z_M as a row vector, and stacked up the results to yield the matrix

then E can be regarded as an incidence matrix³⁰ between the atoms and the quadrature nodes. As shown in Fig. 13, this is similar to the graph-based abstraction as seen in the Coulomb matrix and the GRAPE kernel. However, in both cases, the vertices in the graph represent atoms, while the edges represent pairwise interatomic interactions. Here, the density-based incidence matrix adopts the opposite pattern and constructs a graph with the quadrature nodes being vertices and atoms being edges. The adjacency matrix in this case is computed as the inner product E^TE,

The weight on the edges, as represented by the elements of the adjacency matrix A, can be interpreted as the total flux as contributed by all paths each bridged by an atom k. We have numerically found that the smallest N eigenvalues (except for the 0 eigenvalue) of the symmetric normalized Laplacian

are invariant under rotation up to a certain noise level, even if the quadrature nodes do not rotate with the atoms. Nonetheless, this detour appears to represent a pure theoretical interest rather than any practical value.

In this paper, we presented the Density-Encoded Canonically Aligned Fingerprint (DECAF) by exploring the idea of using smoothed density fields to represent and compare atomistic neighborhoods. One of the key enabling techniques in DECAF is a kernel minisum algorithm, which allows the unambiguous identification of a canonically aligned coordinate frame that can be used for rotationally invariant projection of the density field as well as any associated vector quantities. We have performed detailed analysis to study the behavior of the fingerprint by changing various parameters, such as resolution, smoothing length, and the choice of weight functions. We demonstrate that the fingerprint algorithm can be used to implement highly accurate regressions of both scalar and vector properties of atomistic systems including energy, force, and dipole moment, and could be a useful building block for constructing data-driven next generation force fields to accelerate molecular mechanics calculations with an accuracy comparable to those driven by quantum mechanical theories and calculators.

This work was supported by the Department of Energy (DOE) Collaboratory on Mathematics for Mesoscopic Modeling of Materials (CM4). This work was also supported by the Army Research Laboratory under Cooperative Agreement No. W911NF-12-2-0023.

As candidates for the weight of integral and density scaling functions (Sec. III B), a class of compact polynomials that satisfy the criteria,³¹ 

1. should be compactly supported,

2. should be strictly positive within some cutoff distance r_c,

3. should decrease monotonically,

4. should be at least twice continuously differentiable,

with a minimal number of non-zero terms are

where s = 1 − r/h is the normalized complementary coordinate within the span h of the kernel, and

is an optional normalization factor to ensure that the integral of the kernel in a 3D ball of radius h is unity. The parameters a and b are the free parameters that can be used to adjust the smoothness and width of the kernel and can take any real numbers satisfying the condition a > b > 2. A series of instantiations of the kernel with various $a$⁠, $b$ values are illustrated in Fig. 14. Note that the kernel W^4,3 is equivalent to the Lucy kernel commonly used in Smoothed Particle Hydrodynamics simulations.³² The kernel can be evaluated very efficiently using only multiplication and addition when both a and b are integers.

In Table IV, we list the nodes and weights of the Laguerre quadrature rules up to N_r = 6, using notations from Eq. (22). In Table V, we list the nodes and weights of the Lebedev quadrature rules up to N_r = 6, using notations from Eq. (21). The Laguerre and Lebedev quadrature nodes can be combined using Eqs. (23)–(26) into composite grids for sampling the atomistic density field.