{'fraction_non_alphanumeric': 0.0437734500625335, 'fraction_numerical': 0.01993925317134179, 'mean_word_length': 4.537396121883656, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Accelerated materials property predictions and design using motif-based fingerprints

Data-driven approaches are particularly useful for computational materials discovery and design as they can be used for rapidly screening over a very large number of materials, thus suggesting lead candidates for further in-depth investigations. A central challenge of such approaches is to develop a numerical representation, often referred to as a fingerprint, of the materials. Inspired by recent developments in chem-informatics, we propose a class of hierarchical motif-based topological fingerprints for materials composed of elements such as C, O, H, N, F, etc., whose coordination preferences are well understood. We show that these fingerprints, when representing either molecules or crystals, may be effectively mapped onto a variety of properties using a similarity-based learning model and hence can be used to predict relevant properties of a material, given that its fingerprint can be defined. Two simple procedures are introduced to demonstrate that the learning model can be inverted to identify the desired fingerprints and then, to reconstruct molecules which possess a set of targeted properties.


I. INTRODUCTION
Data-driven approaches towards materials design and discovery are rapidly increasing in popularity, demand and potency. [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15] This emerging trend is fueled by the availability and emergence of large materials databases, [16][17][18] as well as our ability to progressively accumulate materials data via high-throughput computations 19,20 and experiments. [16][17][18] Data-driven strategies aimed at rapid property predictions, and ultimately to rational or informed materials design, rely on exploiting the information content of past data, and using such information within heuristic or statistical interpolative learning models to provide estimates of properties of a new material. This approach is entirely analogous to similar pursuits undertaken within chem-and bio-informatics wherein lead candidates worthy of further in-depth investigations are identified rapidly in a first-level of screening. 4,5,14 Data-driven property prediction strategies have two steps. The first involves representing materials numerically via descriptors, attribute vectors, or fingerprints. In the second step, using available "training" data sets, a mapping is established between the numerical representation of materials and their properties, thus leading to a prediction model. Subsequently, the properties of a new material are estimated using this model after reducing the material to its numerical representation.
One of the central challenges in this whole process is deciding on an appropriate and acceptable numerical representation of materials. The specific choice of this representation is entirely application dependent, and can range from high level descriptors (e.g., d-band center, atomic electronegativities) 21,22 to topological features (e.g., substructural motifs) 20,23,24 to microscopic fingerprints that may capture chemical and configurational degrees of freedom (e.g., coulomb matrix, symmetry functions). [25][26][27][28] Regardless of the specific choice, the representations are expected to satisfy certain basic requirements. These include invariance of the representation with respect to transformations of the material such as translation, rotation, and permutation of like elements. Moreover, it is desired that the representation be intuitive, elegant and physically and chemically meaningful.
In this contribution, inspired by developments in cheminformatics, 14,15 we propose a class of hierarchical motifbased topological fingerprints. This choice, in which the motifs are molecular fragments of varying sizes, is particularly suited to representing molecules and solids composed of elements such as H, C, N, O, F, etc., whose coordination preferences are well understood. Large datasets of molecules and solids are considered, and it is shown that the fingerprints may be effectively mapped to a variety of properties using a similarity based learning algorithm. Moreover, it is demonstrated that the learning model may be inverted to identify fingerprints, and subsequently, to reconstruct actual molecules that possess a desired set of target properties.

II. DATASETS
In the present work, we restrict ourselves to systems composed of C, O and H. We used two datasets, one for molecules and one for crystals, to demonstrate the applicability of the proposed fingerprints. Of these two datasets, the former was taken from Ref. 19 while the latter was prepared by us.

A. Molecule dataset
A dataset of more than 134,000 small molecules made up of C, O, H, N, and F was reported in Ref. 19. This reliable dataset, which contains the optimized geometries, and energetic, electronic, and thermodynamic properties calculated using the B3LYP hybrid exchange-correlation (XC) functional and the 6-31G(2df,p) basis set with the Gaussian 09 software, sets up the stage for many interesting data-mining works. 29,30 A subset of this dataset, containing 45,708 molecules composed of C, O, and H was used in this work. Five properties were considered, including the atomization energy E at , the energy gap E HL between highest occupied and lowest unoccupied molecular orbitals (HOMO-LUMO gap), the isotropic polarizability α, the heat capacity C v , and the zero-point vibration energy E ZP . The obtained structures were optimized by firstprinciples calculations within the DFT formalism as implemented in Vienna Ab initio Simulation Package (vasp), [35][36][37][38] utilizing the semi-local rPW86 XC functional 39 and a plane wave energy cutoff of 400 eV. A Monkhorst-Pack k-point mesh 40 with the spacing of no more than 0.15Å −1 in the reciprocal space were used for sampling the Brillouin zone, while the van der Waals interactions were estimated with the non-local density functional vdW-DF2. 41 Convergence was assumed when the atomic forces exerting on the atomic sites are smaller than 0.01 eV/Å. The entire crystals dataset, which includes the optimized structures, the atomization energies E at , the band gaps E g , and the electronic and ionic parts of the dielectric constants, elec and ion , can be found in the Supplemental Material. 42

III. FINGERPRINTS
A hierarchy of equilibrium structure fingerprints of the same family with increasing levels of sophistication are proposed here. The construction of fingerprints was guided by two simple chemical concepts, i.e., chemical bonds and coordination number. The former intuitively characterizes the short-range interatomic interactions 43 while the latter is the number of bonds involving a given atom. In major classes of materials composed of light elements like C, H, O, N, and F, these concepts are welldefined. In particular, the length of a given bond involving these elements falls in a narrow range (see Refs. 44 and 45 for a comprehensive bond length statistics). For instance, the equilibrium length of a single bond between two C atoms is 1.50Å, the length of a double bond between two C atoms is 1.45Å, and the length of a double bond between a C atom and an O atom is 1.20Å. 44,45 The coordination number is also well-defined, i.e., for a C atom, it can only be 2, 3, or 4 while each O atom can generally bond with 1 or 2 other atoms. Therefore, atoms in a structure can be unambiguously classified (or labeled) by Ai where A is the type of the element (A ∈ {C, O, H}) and i is its coordination number. Likewise, bonds can be specified by the types of its two ends, e.g., Ai-Bj. For the datasets of C, O, and H, the six possible atom types are C2, C3, C4, O1, O2, and H1 while there are sixteen chemically permissible types of bonds, namely C2-C2, C2-C3, C2-C4, C2-O1, C2-O2, C2-H1, C3-C3, C3-C4, C3-O1, C3-O2, C3-H1, C4-C4, C4-O2, C4-H1, O2-O2, and O2-H1. Except C2-O1, C2-O2, and O2-O2, thirteen of them are present in our molecules and crystals datasets. The atom and bond types belong to a family of related structural building units (subsequently described) that can be used to numerically represent the materials structures and hence, are used to define the fingerprints. In particular, the i th −order fingerprint f (i) is defined in terms of its components as Here, n (i) κ is the number of building units (or fragments or motifs) of type κ and N at is the number of atoms either in the molecule or in the unit cell of a crystal. Four types of fingerprints, namely f (0) , f (1) , f (2) , and f (3) , are discussed in the following subsections.
A. 0 th -order fingerprint, f (0) The simplest (0 th -order) fingerprint f (0) represents the fractions of all the element types A existing in the structures, i.e., κ ≡ A. Therefore, in the defini- κ≡A is the number of atoms of element A. This fingerprint is a three-dimensional vector whose components satisfy a simple normalization condi- Next in the hierarchy is the case κ ≡ Ai in which n (1) κ≡Ai is the number of A atoms which are i−fold coordinated. f (1) is a 6-dimensional vector, satisfying several constraints established from the definition or from the chemistry. The first one is the normalization condition, given as Within the two datasets, all the C2 atoms should be grouped by pairs, forming triple C≡C bonds. Therefore, the number of C2 atoms, which is N at × f C2 , must be an even integer. Moreover, since each C3 atom only make a double bond with either an O1 atom or another C3 atom, one must have f O1 is an even number. By examining the connectivity of a structure, another constraint reads where N is the number of closed loops of bonds and d is a structure-dependent parameter. For molecules and crystals composed of isolated substructures (or molecules), d = 0 while for crystals composed of connected substructures, d > 0. The derivation of this constraint is given in Appendix A. The last constraint of f (1) is written in the form of a recursion relation, i.e., C. 2 nd -order fingerprint, f (2) Both f (0) and f (1) are local, representing the density of the atom types of a material. The equilibrium interatomic distance is somehow captured by the 2 nd -order fingerprint f (2) where all the possible bonds are counted. f (2) is a 13-dimensional vector whose components, f Ai-Bj , represent the normalized number n (2) Ai-Bj of the Ai-Bj bonds in the structure. From f (2) , f (1) can readily be determined by a recursion relation where δ Ai,Bj is used to remove the double counting when Ai ≡ Bj [see Appendix B for the derivation of (5)]. Through this recursion relation, all the constraints that f (1) obeys are applicable for f (2) . We note that f (2) was discussed in several previous works, e.g., in Refs. 25, 46, and 47 under the name of "bond counting". This fingerprint can also be regarded as a generalization of "doubles", the fingerprint defined in Ref. 20 for the chain models of polymers.
In the 3 rd -order fingerprint f (3) , the number of twobond catenation is represented, i.e., κ ≡ Ai-Bj-Ck. In particular, the definition (1) for f (3) κ≡Ai-Bj-Ck involves n Ai-Bj-Ck , which is the number of Ai-Bj-Ck sequences, or equivalently, the catenation of two bonds Ai-Bj and Bj-Ck. Considering compounds of C, O, and H, there are 125 possible distinct catenation of two bonds Ai-Bj Similar to f (2) , f (3) can be viewed as a generalization of "triples", the fingerprint examined in Ref. 20.

IV. PROPERTY PREDICTION MODEL
A learning model is critical in order to map the fingerprints to properties. In this work, we chose Gaussian kernel ridge regression (KRR), 5,48,49 the technique which has successfully been used in material properties predictions 20,25,[28][29][30] Within this model, the input fingerprints are transformed into higher-dimensional space whereby a linear relation between the transformed fingerprints and the associated properties can be established. This mapping involves the distances between fingerprints and can be regarded as a similarity-based prediction model, i.e., similar properties may be predicted for materials with similar fingerprints.

FIG. 2. (Color online)
Learning curves corresponding to Eat, EZP, α, Cv, and EHL. For each model, (2) , and f (3) are used to represent the molecules. Calculated data is given by symbols while curves are the guide for the eyes.
where ν runs over all the fingerprints in the training dataset. Here, d µν is the distance between fingerprints µ and ν, defined as the Euclidean metric d µν = The Gaussian width parameter σ and the regression coefficients α ν are determined within the training phase whence a regularized objective function is minimized. 5,48,49 During this phase, σ and the regularization parameter are determined by k-fold cross validation on the training set (k = 5 in this work). Within this method, the training dataset is split into k bins, any of the bins is considered to be a new test dataset while the remaining k − 1 bins form a new training datatest. This procedure is repeated for each of the k bins and for every value of σ and λ on a preselected logarithmic-scale grid.
The optimal values of σ and λ, i.e., those leading to the minimum k-fold cross-validation (mean absolute) error, are used to compute α ν of the entire dataset.

A. Molecules dataset
The four fingerprints considered, namely f (0) , f (1) , f (2) , and f (3) , were used to represent the molecules dataset. To mimic the learning and prediction processes, the dataset was randomly partitioned into a training dataset and a test dataset. The KRR model was then trained on the training dataset using five-fold cross validation before predictions were made on the test dataset. We show in Fig. 2 the learning curves of E at , E ZP , α, C v , and E HL , plotting the training and test errors against the number of molecules in the training dataset (data reported in this figure was averaged over 30 independent runs). In addition, predictions for the test dataset of 44,708 molecules after training the KRR model on a dataset of 1,000 molecules are shown in Fig. 3. As discussed in detail below, both Fig. 2 and Fig. 3 indicate that all of these properties can be very well predicted by using either f (2) or f (3) , provided that the KRR model is trained on a training dataset of 200 or more data points.
The general tendency, as revealed by Fig. 2, is that higher-order fingerprints offer more accurate predictions. The 0 th -order fingerprint f (0) can be used to roughly estimate energy-related quantities, i.e., E at and E ZP while it can not be used for others. For instance, E HL can not be predicted with f (0) because this fingerprint is totally local in nature, encoding no information at any finite range. Consequently, the finite conjugation length, known to signal the energy gap reduction in complex (conjugated) systems (see, for example Ref. 50), is not captured by f (0) . Fingerprints of higher orders, e.g., f (1) , f (2) and f (3) , contain some information at increasing ranges, allowing for systematically better predicting E HL . These fingerprints also work sufficiently well in predicting E at and E ZP . With f (1) , the averaged error in predicting E at is 25 meV/atom while this error is reduced to 20 meV/atom and 18 meV/atom if f (2) and f (3) , respectively, are used. The very good power of f (2) in predicting E at reproduces the similar conclusions drawn for the "bond counting" fingerprint by Ref. 47. This behavior is understandable because the dissociation energy of chemical bonds in organic molecules and crystals, which dominates the stability of these systems, are well-defined 46 in the same fashion with the bond length as previously discussed. Interestingly, this predictive power can significantly be improved if more advanced fingerprints, i.e., those can capture the small perturbations of interatomic distances like Coulomb matrix, are used. 29,30 Compared to f (1) and f (2) , f (3) is significantly better in predicting C v . The considerable improvement in the predictions of α when f (2) is used instead of f (1) may indicate the key contribution from polar bonds to the high-value regime of α.

B. Crystals dataset
We performed similar predictions for the dataset of 215 crystals containing C, O, and H. Using the KRR model coupled with f (0) , f (1) , f (2) , and f (3) , five properties of these crystals, including the atomization energies E at , the band gap E g , the electronic dielectric constant elec , the ionic dielectric constant ion , and the total dielectric constant tot = elec + ion , were predicted. We show in Fig. 4 the learning curves, representing the errors of the predictions using these fingerprints, averaged over 100 independent runs. In Fig. 5, the predictions for the five properties are given, using the KRR model trained on a random training set of 150 data points.
Clearly, the tendency of the prediction performances on the crystals dataset is similar to those of the molecules dataset, i.e., high accuracies are obtained with fingerprints of higher orders, and properties which are governed by long-ranged information, e.g., band gap E g , can only be predicted with high-order fingerprints. For the atomization energy E at , predictions with f (0) and f (1) leads to quite high averaged errors, which reduced to 18 meV/atom and 15 meV/atom when f (2) and f (3) , respectively, were used. Overall, all the five examined properties can be predicted well when high-order fingerprints are used to represent the crystals. For instance, by employing f (3) , the averaged error in predicting E g is 0.45 eV while the electronic dielectric constant elec and the ionic dielectric constant ion can be predicted with an averaged error of 0.1 − 0.2.

VI. UTILITIES OF THE FINGERPRINTS
The demonstrated predictive power of the KRR model, which uses f (i) to represent materials structures, inspires the idea of using this model to rationally optimize materials for a targeted property P opt , the concept often referred to as "inverse design". [51][52][53][54] In fact, a large number of success stories along this direction have been reported in the past, using various approaches, e.g., iteratively optimizing the properties of a given compound or on-the-fly screening when searching for stable structures. [55][56][57][58][59][60][61][62][63][64][65][66][67] Here, our idea is that starting from a trained KRR model, fingerprints which correspond to the desired properties can be predicted. Then, molecular structures will be reconstructed from the predicted fingerprints. Finally, the targeted properties will be verified by DFT calculations at the same level with those used for the training dataset.
The greatest challenge of this procedure is to ensure that the predicted fingerprint is physically and chemically meaningful, i.e., at least one material structure can be reconstructed from it. 68,69 Therefore, one must mathematically define the subspace of the meaningful fingerprints, and then limit the search for desired fingerprints within this subspace. We present two approaches which can be used for designing molecules (the work of designing crystals is not considered here).

A. Design via enumeration
The central idea of this approach is that the components of a given fingerprint can be enumerated in a given way so that it is meaningful. We used f (2) for a demonstration because predictions using this fingerprint are good while its dimensionality is not too high like f (3) . We first implemented the applicable rules involving bonds and coordination numbers by defining five "backbone" blocks. They include C4, C=C (a pair of C3 atoms with a double bond), C≡C (a pair of C2 atoms with a triple bond), C=O (one C3 and one O1 atom linked by a double bond), and O2. By definition, all of the dangling bonds starting from these blocks are single, thus any of them can be connected to others without any constraint. Then, given a set of backbone blocks, all the possible arrangements can be scanned, keeping track of the connectivity to eliminate some dangling bonds, and saturating the remaining dangling bonds by either H1 or OH, referred to as "ending" blocks. From the obtained arrangements, f (2) can be unambiguously deter-mined and their properties were predicted. Those with targeted properties were singled out to rebuild molecular structures for validating calculations. We show in Fig.  6 two optimized molecules constructed from two of the predicted fingerprints, labeled by A and B, accompanied by the predicted and calculated E HL and α. The results given in Fig. 6 indicate that the desired molecules are indeed obtained.

B. Design via inversion
Different from the enumeration approach, this procedure aims to directly determine the fingerprints, starting from desired properties. This goal can be achieved by optimizing an objective function, aiming towards the desired properties while applying the constraints that ensure the fingerprints considered are meaningful. Because the reconstruction step requires a simple enough fingerprint, f (1) was selected for this approach. Among the constraints established for f (1) , (2) and (3) are explicitly imposed in the objective function defined below Here, λ 1 , and λ 2 are the Lagrange multipliers associated with the constraints while P is the property (or properties) of the trial fingerprint f (1) predicted by the trained KRR model. In practice, we evaluated P by averaging many predictions, each of them was given by the KRR model trained on a randomly selected training dataset of 1,000 data points. All the terms in (8) are given in the quadratic form to smoothen G. Generally, the problem of minimizing G[f (1) , λ 1 , λ 2 ] (performed with simulated annealing 70 in this work) returns many solutions F (1) . For each of them, N at was determined by minimizing another objective function D[F] defined as where nint(x) returns the closest integer to x. Once N at is determined, a post-screening step is performed to consider the possibility of N > 0 and to single out the fingerprints so that N at F (1) C2 and N at F are positive even numbers. Such fingerprints are meaningful, i.e., molecules can be built up from any of them.
We demonstrate this procedure by optimizing two properties simultaneously, i.e., E HL and α. We note that these properties seem to be competing, as shown in Fig. 7 where an asymptotic limit of the form α ∼ 1/E HL can be seen (similar limit between two related properties of crystals, namely elec and E g was documented earlier in Ref. 71). An examination of Fig. 3 reveals that the prediction of α using f (1) is fairly good in the region of α < 0.8 eV, all the predictions for E HL by the KRR model are below 9 eV. The reason is given in Fig. 3 which clearly implies that when f (1) is coupled with the KRR model, high values of E HL (8 ≤ E HL ≤ 10 eV) are generally underestimated by roughly 1 eV. Three of the predicted fingerprints, labeled by C, D, and E, were selected for rebuilding new molecules. From either C or E, only one molecule can be constructed while many different molecules correspond to D. All of the molecules reconstructed from C, D, and E were optimized and then their α and E HL were calculated with Gaussian 09, 72 using the 6-31G(2df,p) basis set and the B3LYP XC functional. 73,74 The results are summarized in Table I

C. Remarks
It is worth noting that the key feature of f (i) which is useable for the described enumeration and inversion design procedures is their discontinuity with respect to slight configurational perturbations. Because all the possible chemical bonds appearing in a molecule comprising C, O, and H are well-defined, it is very likely that the optimization step performed on the reconstructed molecules preserves the predicted fingerprint. Moreover, the efficiency of the designing approaches depends on several factors, including the prediction accuracy of the fingerprints used. Although predictions by using highorder fingerprints are systematically better, the complexity generated by their high dimensionality is significant. Comparing to the procedure described above, that utilizing f (2) or f (3) needs roughly 10 and 100 more constraints for ensuring the considered fingerprints are meaningful.
If the dimensionality of f (2) can considerably be reduced, it may then be used for the inversion approach.

VII. CONCLUSIONS
To summarize, we have systematically studied a family of motif-based topological fingerprints which can numerically represent major classes of molecules and crystals. By using a similarity based learning algorithm, these fingerprints can be mapped onto various properties of molecules and crystals, significantly accelerating their properties prediction. A major advantage of these fingerprints is clearly demonstrated via two procedures for designing molecules, one by enumeration and the other by inversion. These procedures rely on the accelerated properties prediction to identify the desired fingerprints, and then to reconstruct molecules that possess one or more targeted properties. We note that although only molecules and crystals comprising C, O, and H are considered in this contribution, our results can straightforwardly be generalized to those containing other light elements whose coordination preferences are well established, e.g., N and F.

ACKNOWLEDGMENTS
The authors thank Venkatesh Botu, Ghanshyam Pilania, and Vinit Sharma for useful discussions and O. Anatole von Lilienfeld for drawing our attention to some important relevant works. The present work was supported by a Multi-University Research Initiative (MURI) grant from the Office of Naval Research, under award number N00014100944. Part of the computational work was done with our sponsored TeraGrid XSEDE allocation. 75 Appendix A: Constraint of f (1) derived from elementary chemical rules Constraint (3) was derived with an assumption that the desired molecular structure is connected, i.e., any pair of atoms are connected by at least one sequence of the allowed chemical bonds. Let us take a molecule in which n Ai is the number of the blocks Ai. Starting from the applicable chemical rules, all the two-fold coordinated carbon atoms are grouped by pairs, forming n C2 /2 units of C ≡ C, each of which is a pair of carbon atoms linked by a triple bond. Next, n O1 one-fold coordinated oxygen atoms must bond with n O1 three-fold coordinated carbon atoms to form n O1 units of C = O. Then, the remaining n C3 − n O1 three-fold coordinated carbon atoms are grouped together by pairs, forming (n C3 −n O1 )/2 units of C = C. Therefore, the set of the blocks Ai now contains n C2 /2 + n O1 + (n C3 − n O1 )/2 + n C4 + n O2 units of C ≡ C, CO, C = C, C4 and O2. Assuming that these units are isolated, the total number of dangling bonds starting from them is 2(n C2 /2)+2n O1 +4[(n C3 −n O1 )/2]+4n C4 + 2n O2 , or simply n C2 + 2n C3 + 4n C4 + 2n O2 . (A1) By joining n C2 /2+n O1 +(n C3 −n O1 )/2+n C4 +n O2 units together, the number of dangling bonds that will be annihilated to form inter-unit bonds is 2[n C2 /2+n O1 +(n C3 − n O1 )/2 + n C4 + n O2 − 1] + 2n where n is the number of loops of bonds, each of which costs extra 2 bonds. Therefore, the number of remaining dangling bonds is n C3 + 2n C4 − n O1 − 2n + 2. (A2) All of these dangling bonds must be saturated by n H1 hydrogen atoms, thus n H1 = n C3 + 2n C4 − n O1 − 2n + 2.
The constraint (3) can then be obtained when we divide Eq. (A3) by N at . This constraint is applicable not only for molecules but also for crystals formed by repeatedly placing an isolated molecule in a periodic grid. If these molecules are not isolated, i.e., they form a network of d dimensions, 2d dangling bonds are used to form the network (assuming that the network are formed only by single bonds). Thus, Eq. A3 is given as In the general case when not only single bonds involve the network formation, the parameter d used in Eq. A4 is not necessarily an integer.