Chapter 1.4. Data mining. I. Machine learning in crystallography

In crystallography (as in all other fields) obtaining and handling experimental data (including data collection, data cleansing and data warehousing) is usually the first, and always a very important, stage. The comprehensiveness of the structural data collection, the structure and quality of the data sets, and the selection of relevant data sets are extremely important for getting reasonable results from the fully automatic procedures of data mining. The application of data mining in crystallography is natural because crystallographic data sets are very large. They are stored in the three main databases: the CSD, the main repository for small-molecule organic and metal-organic crystal structures (more than 1 160 000 entries); the ICSD (completely identified inorganic crystal structures; more than 185 000 entries); and the PDB (biological macromolecular structures; more than 187 000 structural hits). The Powder Diffraction File from the International Centre for Diffraction Data, with over 1 076 000 powder diagrams, should be mentioned too. For discussions of these databases, we refer the reader to numerous reviews and publications, for example Buchsbaum et al. (2010), Allen (2002) and Berman (2008). Recent developments in CSD tools have been summarized by Groom et al. (2016).

The collections of experimental data serve as input information in the data-mining procedure. In the data-mining literature one can also find the terms `input', `input object', `input vector', `input variables' and so on. In the following we will use the term input vector, . The input vector in data mining can contain the values of the various descriptors¹ for an individual object, which can be given their own weights w (in the following the terms weight vector and weight matrix are used), reflecting the importance of a given descriptor for a given input vector. The input vector can be presented as a list of values or as a row in a table (matrix). In a more general way, the `descriptors' are the template for defining individual objects. In crystallography the individual objects are obviously the crystals, and the most widespread template for describing individual crystals and their structures is the Crystallographic Information Format (CIF) [see International Tables for Crystallography, Volume G (Hall & McMahon, 2006)]. The CIF template contains fields for the unit-cell parameters, crystal symmetry, atomic coordinates, molecular formula, crystal density, data on the X-ray experiment, temperature, pressure and so on. Which of these descriptors will be used in the input vector depends on the target property or pattern of interest. To complete the input vector, data filtering is normally applied. Here it is worthwhile to note that in the following discussions we will distinguish between the filtering and the cleansing of data. While the filtering serves to extract subsets from the whole data set according to some features, for example the experimental conditions, cleansing is the process of detecting and correcting (or removing) incomplete, incorrect, inaccurate or irrelevant parts of the data.

Data analysis and model development are consecutive, and strongly correlated, steps that use the algorithms of machine learning. The choice of appropriate instruments and approaches (or their combination) depends on the current task. The global goals of `pattern recognition' (description) and `pattern prediction', mentioned above, are normally achieved in crystallography within the following most important data-mining tasks:

Clustering – This is the task of discovering groups in the data that were not predefined. The members of each group are `similar' in some way. Cluster analysis, when integrated with the proper similarity function, can serve as a powerful tool for crystal structure prediction, crystal structure determination and rapid screening of databases. The applications of cluster analysis in crystallography will be discussed in Section 3.

Anomaly detection – This does the opposite to clustering, by identifying data or observations significantly different from the norm. The ability to detect such anomalies (or outliers) is critical in materials science because it can help to identify a new class of materials with an unusual property, or to anticipate potential harmful effects which would normally only appear through a retrospective analysis. Anomaly detection can be used to isolate faulty data in databases and clean initial data sets. It can identify unusual data records that might be interesting and reward further investigation. Some examples of uses of anomaly detection in crystallography will be given in Section 5.

Classification – This is the task of finding a function that arranges the data into predefined groups that may have well behaved correlations and can form the basis of physically and statistically based models. Classification allows estimation of physical properties as well as the derivation of model parameters (for example, parameters of force fields). Classification predictive modelling approximates a mapping function from input vectors to discrete output vectors (for more details see Section 4).

Regression – Regression attempts to find a function that maps the input vectors onto predefined target values. These target values are the measurement of a physical property, which should be predicted by the parameterized model. Regression predictive modelling maps the obtained function from input vectors (variables) to a continuous output vector. The success of the regression models also helps to refine the usefulness and relevance of the input vectors. Some results obtained with regression will be given in Section 4.

Summarization (dimensionality reduction) – This provides a more compact representation of the data set, including visualization and report generation (Section 6).

From a mathematical point of view, the tasks and the algorithms to manage the tasks can be classified as unsuper­vised or supervised methods. Cluster analysis is an unsuper­vised method, whereas classification and regression are supervised methods. Anomaly detection and summarization can be either (or both). For example, summarization (also known as dimensionality reduction) is the process of reducing the number of random descriptors (also referred to in the literature as `variables', `attributes' or `features') needed to describe the variations in the data set. In some cases, data analysis such as regression or classification can be done in the dimensionally-reduced space more accurately than in the original space. It can be unsupervised (e.g. principal component analysis, PCA, and singular value decomposition, SVD), or supervised (e.g. linear discriminant analysis, LDA). The algorithms most often used in crystallography are recursive fitting, decision tree, support vector machines (SVMs), singular value decomposition (SVD), neural networks (NNs) and genetic algorithms (GAs). Some of them have specific applications, others are broader in application and can be used in both unsupervised and supervised tasks. For example, recursive fitting was used in Hofmann et al. (2007) for deriving a `reactive' potential for water. SVMs for classification tasks have been used and described in detail in a series of works (Hofmann & Apostolakis, 2003; Apostolakis et al., 2001; Hofmann, 2010). The application of SVD will be described in Sections 4.2 [see also Hofmann (2002)] and 4.3.4.

Neural networks (NNs) were first described by McCulloch & Pitts (1943) and achieved enormous popularity in the 1990s (Craven & Shavlik, 1997; Boser et al., 1992), especially in protein crystallography [see, for example, Riis & Krogh (1996)]. The neural-network algorithm was attractive owing to its high flexibility and ability to fit any type of function. However, the choice of network architecture is not always trivial and therefore the method often produces incomprehensible models and requires long training times in classification tasks (Maimon, 2010; Zhang, 2010; Nikam, 2015). These days, these types of problems are addressed in a significantly simpler way by SVMs (Apostolakis, 2010). Even so, both algorithms (SVMs and NNs) are used in data mining to manage classification and regression tasks.

It is important to note that SVD allows for the physical interpretation of results (weight vectors, mathematical models and parameters, Section 4.3.4), while for NNs this is not possible owing to the large number of parameters in the concatenated matrix (Fig. 2).

Figure 2 Generalized NN flowchart for predicting the structure of ice. The sum formula is multiplied by the weight matrix to propagate from one layer to the next layer. The multiplication in the last layer produces the desired crystal-structure properties as the output vector.

The final step of knowledge discovery from data is to verify that the patterns produced by the data-mining algorithms occur in the wider data set. Not all patterns found by the data-mining algorithms are necessarily valid. It is common for the algorithms to find patterns in the training set that are not present in the general data set. This is called `overfitting', and often happens if a small number of data points are used. To overcome this, the evaluation uses a test set of data that the data-mining algorithm was not trained on. The learned patterns are applied to this test set and the resulting output is compared with the desired output. A number of statistical methods may be used to evaluate the algorithm, such as receiver operating characteristic curves. If the learned patterns do not meet the desired standards, then it is necessary to re-evaluate and change the preprocessing and data mining. If the learned patterns do meet the desired standards, then the final step is to interpret the learned patterns and turn them into knowledge.

If databases also contain target values, or if a measurement for the desired target property exists for at least some input vectors, the supervised learning algorithms `learn' how to predict these values with new records when the output is unknown. More generally, they can also consider several properties of interest. In this context `the target property (pattern)' might be more appropriately defined as the output vector or output variables. In the data-mining literature the values being predicted are also called the `dependent variables', `output variables', `target variables' or `outcome variables'.

Data clustering according to similarity is a common technique for statistical data analysis. It is used in many fields, including machine learning, pattern recognition, image analysis and bioinformatics (Duda & Hart, 1973). Clustering is normally performed using one of two types of clustering algorithms; hierarchical or partitional (Fig. 3). Hierarchical algorithms find successive clusters using previously established clusters, whereas partitional algorithms determine all clusters at once. Partitional algorithms divide a set, , into non-overlapping `parts' (or `blocks' or `cells') that cover all of , Fig. 3(a). More formally, these cells are both collectively exhaustive and mutually exclusive with respect to the set being partitioned. Hierarchical algorithms can be agglomerative (`bottom-up') or divisive (`top-down'), Fig. 3(b). Agglomerative algorithms begin with each element as a separate cluster and merge them into successively larger clusters. Divisive algorithms begin with the whole set and proceed to divide it into successively smaller clusters. Two-way clustering, co-clustering or bi-clustering, which should be mentioned here, are the names for clustering where not only the objects but the features of the objects are clustered (i.e., if the data are represented in a data matrix, the rows and columns are clustered simultaneously).

Figure 3 (a) Partitional and (b) hierarchical clustering.

The most popular algorithm for clustering in biology, pharmacy and crystallography is hierarchical clustering, with the traditional representation of the hierarchy as a tree, or a dendrogram, Fig. 3(b), with individual elements at one end and a single cluster containing every element at the other. Familiar examples of dendrograms are family trees and the evolutionary tree of Darwin. Cutting the tree at a given height will give a clustering at a selected precision, e.g. in family trees different levels of relationship can be selected. Cutting at zero gives identities.

The common algorithms for agglomerative hierarchical clustering are single linkage, complete linkage, and average or centroid linkage. These methods differ in how the distance (similarity) between each cluster is measured. Clustering is started by assigning each structure of the set its own cluster. For example, if we have set of 50 crystal structures, we will define 50 different clusters. Clusters are then successively merged according to the calculated distance (similarity) between elements of the clusters. The distance between two clusters with elements a_i and with elements b_j (Fig. 4) is defined by minimum, maximum or average values of a similarity s_ij (distance) (Apostolakis, 2010) between elements of different clusters as follows.

Figure 4 Hierarchical clustering according to similarity.

Single linkage. In single-linkage clustering, the distance between two clusters and is equal to the minimal (shortest) distance between elements a_i and b_j in each cluster: It is therefore very efficient and simple to understand; it tends, however, to produce large elongated clusters that are not necessarily internally compact. The reason for this is that the merging criterion does not take the overall structure of the new cluster into account.

Complete linkage. In complete-linkage hierarchical clustering, the distance (similarity) between two clusters and is defined as the maximal (longest) distance between elements a_i and b_j in each cluster: Complete linkage tends to create many small clusters at intermediate levels of the cluster tree. While complete-linkage clustering is not able to identify clusters of complex (e.g. elongated) structure, even when they clearly exist in the data, it has often been observed in practice to produce very useful hierarchies.

Average or centroid linkage. In average-linkage hierarchical clustering, the distance between two clusters and is defined as the average distance between each element a_i in one cluster to every element b_j in the other cluster: where and are the number of elements in the clusters. While average-linkage clustering can overcome the problems of the other two methods, it requires more computation. In some versions the average distance to the centroid is taken, which then again corresponds to a minimization of the squared error, in the sense that the two clusters are merged that lead to the lowest possible (average) deviation from the cluster representative (the centroid).

Each agglomeration occurs at a greater distance between clusters than the previous agglomeration, and one can decide to stop clustering either when the clusters are too far apart to be merged (distance criterion) or when there are a sufficiently small number of clusters (number criterion). A graphical representation of the results of hierarchical cluster analysis is the dendrogram, Fig. 4. The branches represent clusters obtained at each step of hierarchical clustering.

We also mention k-means clustering, a simple partitional algorithm which provides good results when the number of required clusters is predefined, or the user has a reasonable idea what it is (Rajan, 2010). The results of k-means clustering are often presented as `structural maps'. The principal components obtained from PCA may be used as indicators for k-means clustering (Ding & He, 2004); see, for example, Balachandran & Rajan (2012) and Rajan (2010). This approach is widely used for crystal structure (and property) prediction in inorganic crystallography.

The calculation of the similarity indexes and subsequent clustering can be used for rapid direct screening of databases to find polymorphic and isostructural crystal structures and manifold entrances. This is because we have to take into account only a very limited number of reflections. The calculation time depends on the choice of the interval θ_min and θ_max. For simulations of powder diagrams the interval from 0° to 60° was used in the work we describe here. The upper value is large enough to include reflections even for very small cells, and there is no limit for the lower value in simulated powder diagrams. In experimental powder diagrams at small angles there is a blind region and the minimum angle is greater than zero.

Specific examples of cluster analysis help us to visualize and to understand the effects of external conditions on structural changes in crystals.

As an example, we show the clustering of benzene polymorphs, which allowed the detection and analysis of structural changes in polymorphs caused by temperature and pressure variations (Hofmann & Kuleshova, 2014). The set of 15 structures was obtained by extracting the experimental crystal structures with atomic coordinates from the CSD. The results of clustering are presented in Fig. 7. The clustering threshold is set at s(i, j) = 0.012, corresponding to the RMSD similarity criterion used for the Cambridge Blind Test (Day et al., 2009; Bardwell et al., 2011; Reilly et al., 2016). The clustering clearly revealed high-pressure polymorphs II (16, 17) and III (03, 04), which exist over a narrow pressure interval. The entries for polymorph I, despite being defined over a broad range of temperatures and pressures, are subdivided into three structural groups reflecting the external pressure and temperature conditions: a sub-cluster consisting of entries 00 (indicated just as BENZEN in Fig. 7), 02, 15 and 18 (defined at quasi-ambient conditions), a single structure 12 (at extreme high pressure) and a cluster comprised of entries 01, 07, 06, 14, 13 and 11 (defined at non-ambient conditions). The crystal densities in the quasi-ambient and non-ambient clusters differ significantly: from 1.02 to 1.12 g cm⁻³ in the first cluster and from 1.14 to 1.26 g cm⁻³ in the second cluster.

Figure 7 Clustering of benzene polymorphs reflects structural changes caused by temperature and pressure variations.

Clustering followed by summarizing of results is a crucial instrument in CSP, where a huge number of generated crystal structures have to be analysed for similarity, and multiple and manifold structures have to be removed. A simple ranking based on an energy criterion does not work here, because, firstly, several hundred crystal structures can be found for one molecule within a very narrow energy range and, secondly, completely different crystal structures can have nearly identical energies.

As an example, we analyse the results of CSP for alanine, an α-amino acid that is used in the biosynthesis of proteins. It contains an α-amino group (which is in the protonated form), an α-carboxylic acid group (which is in the deprotonated COO⁻ form) and a side-chain methyl group, classifying it as a nonpolar (at physiological pH) aliphatic amino acid. The L-isomer (left-handed version) of alanine is one of the 20 amino acids encoded by the human genetic code. The right-handed form, D-alanine, occurs in bacterial cell walls and in some peptide antibiotics. CSP was performed for this molecule with FlexCryst (Hofmann et al., 2016). The CSP algorithm used by this program is described in detail in Hofmann (2010). It is implemented with a particular energy function (Hofmann & Apostolakis, 2003) derived by data mining. This force field is able to reproduce crystal structures with hydrogen bonds (Hofmann et al., 2004) and special sorts of atoms, including ions (see Sections 4.3.4 and 4.4). Because the reference molecule is optically active, the calculations were performed with the most frequently occurring acentric space groups P2₁ and P2₁2₁2₁. After 10 minutes, 200 structures had been generated with a coverage (completeness of crystal structure generation) of 76% (Hofmann et al., 2009). In the following steps the generated crystal structures then needed to be clustered, manifold structures needed to be removed and the resulting structures needed to be compared with the experimental one.

For alanine, additional problems arise: the CSD contains in total more than 90 entries for crystal structures of alanine. Among them are 48 entries for L-alanine (eight defined at ambient conditions of temperature and pressure), 21 entries for D-alanine (one at ambient conditions) and 23 entries for DL-alanine (three at ambient conditions). So, firstly, the reference experimental structure had to be chosen. In Fig. 8 we show the clustering of crystal structures that were only determined at ambient conditions. At the threshold s(i, j) = 0.012, two clusters were defined: one cluster with the DL form of alanine (DLALNI) and the other with optically active L and D forms (LALNIN and ALUCAL). Since the similarity works in reciprocal space, the L and D forms fall into one cluster. After the summarization, two structures were retained: ALLUCAL04 (space group P2₁2₁2₁) and DLALNI01 (space group Pna2₁). Both of them were used as reference experimental structures of alanine.

Figure 8 Clustering of experimental structures of alanine from the CSD. Two structures (marked with a cross) are retained as representative after summarization.

The clustering of generated (predicted) crystal structures is shown in Fig. 9. For clarity of visualization, we present only the 20 top-ranked crystal structures with the lowest values of crystal energy. At the threshold [s(i, j) = 0.012], four clusters were retained; from each cluster only the structure of lowest energy is retained as representative after summarization.

Figure 9 Clustering of predicted crystal structures of L-alanine (the 20 top-ranked structures are presented). From each cluster only the structure of lowest rank (marked with a cross) is retained after summarization.

Finally, the experimental and predicted crystal structures were compared (Fig. 10). At the standard threshold of s(i, j) = 0.012, the correctly predicted crystal structure was revealed as that of rank 15.

Figure 10 Comparison of the retained predicted structures with the experimental structures revealed the correct crystal structure to be of rank 15.

Comparing crystal structures within the molecular coordination spheres consisting of the 15 closest neighbours (the criterion in RMSD), performed with Mercury (CCDC, 2016b), gave the same results. In Fig. 11 we show the overlay of the predicted crystal structure of rank 15 and the experimental crystal structure ALUCAL04.

Figure 11 Overlay of the experimental (grey) and predicted (green) structures of alanine, r.m.s.d. = 0.272.

The method presented here was used for the evaluation of the results in the fourth blind test organized by the Cambridge Crystallographic Data Centre (Day et al., 2009), along with the traditionally used estimations of r.m.s.d. values using RMSD (Chisholm & Motherwell, 2005) and of energy functions (van Eijck, 2005). The results obtained with all three methods were identical. However, the cluster method using similarity introduced by Hofmann & Kuleshova (2005), working in reciprocal space, reduced the calculation times by orders of magnitude.

The possibility of easily finding isostructural entries in databases hints at the use of clustering for crystal structure determination from powder diagrams. In fact, if the available single crystals are too low quality for single-crystal structure determination and only an experimental powder diagram can be obtained, it is worth performing a search of databases for entries that might be isostructural. If an isostructural entry is found in a database, it can be used as a starting point for finding the crystal structure corresponding to the experimental powder diagram. In this way, the crystal structures of Pigment Red 3 (4-methyl-2-triphenylarsine-2-naphthol, CSD refcode MNIPZN) and Pigment Orange 5 [1-((2,4-dinitrophenyl)azo)-2-naphthol, CSD refcode CICCUN] (Schmidt et al., 2007) were solved. In a similar way, the crystal structure Pigment Red 170 was solved: first the crystal structure of the methyl derivative of Pigment Red 170 was found, and then the crystal structure of Pigment Red 170 was found to be isostructural to it (Schmidt et al., 2006).

Another example is crystal structure determination from unindexed powder diagrams. The algorithm for this works as follows: first the crystal structures are generated (using CSP), then their powder diagrams are simulated and compared (clustered) with the available experimental powder diagram. Afterwards, the structures most similar to the experimental crystal structure are minimized according to the evaluation function (in this case it is the similarity between the experimental and simulated powder diagrams). In principle, crystal structure determination from powder X-ray data can be considered as a combination of CSP and clustering of a virtual data set. In the case of indexed powder diagrams, the X-ray diffraction patterns can be included as additional information during crystal structure generation as a boundary condition on the lattice parameters.

The advantages of the integral similarity index for crystal structure determination are (Hofmann & Kuleshova, 2005): similar powder diagrams can be recognized, even if peaks are shifted; the similarity index is calculated quickly because the evaluation works pointwise; and during refinement, interchanged peaks can be reordered correctly. The most widespread similarity index used for comparing powder diagrams is the mean-square deviation between two diagrams (R value, as used in Rietveld refinement). Unfortunately, this metric is very sensitive to even very small shifts of the peaks, and does not allow minimization to be performed if the unit cells have slightly different parameters. The similarity indices mentioned above (Karfunkel et al., 1993; Hofmann & Kuleshova, 2005) are tolerant of small cell differences. However, the index proposed by Karfunkel et al. (1993) arbitrarily provides either positive or negative values of similarity and therefore cannot be used for refinement.

As an example of a crystal structure determination, we show the case of Pigment Red 181 (Hofmann & Kuleshova, 2005). The comparison of simulated powder diagrams of predicted crystal structures with the experimental powder diagram and subsequent structure refinement indicated that the candidate of rank 12 was the most promising crystal structure. The initial value of the similarity index of this structure was 0.036, and this decreased (improved) during refinement to 0.003.

In Fig. 12, one can observe shifted peaks at 10 and 20° 2θ for the predicted crystal structure. The shift of these peaks (the 010 and 020 peaks) is caused by the predicted cell parameter being 0.4 Å too short. During refinement, this inaccuracy was corrected and the cell parameter was elongated from 9.21 to 9.63 Å. Another problem for refinement would be caused by the two peaks at 25° and 27°. In the powder diagram for the predicted crystal structure, these peaks overlap and are located at 27°. The algorithm correctly recognizes that this peak has to be split and separates it. As for the case of interchanged peaks, other similarity indices would have a local minimum for this situation, making further refinement impossible. The packing of the molecules that was determined is shown in Fig. 13.

Figure 12 Structure determination of Pigment Red 181. From the predicted polymorphs, the structure giving a simulated powder diagram closest to the experimental powder diagram is selected. After refinement, the powder diagrams are nearly identical. Reprinted from Hofmann & Kuleshova (2005).

Figure 13 Molecular packing of Pigment Red 181. Reprinted from Hofmann & Kuleshova (2005).

The goal of supervised machine learning is the parameterization of a model function, f, based on labelled training data P^exp, to produce (simulate) the desired property P^sim as well as possible. Ideally P^exp = P^sim. The training data consist of a set of training examples. In supervised learning, each example k is a pair consisting of an input object (typically an input vector ) and an output value, . The labelled input objects are different descriptors. In our cases they contain the crystallographic, chemical and some physical properties of the compounds, for example the sum formula, all atom–atom distances, the atom weights or the sublimation heat. The desired output value might be the crystal density, the crystal structure, the solubility, the sublimation heat, the chemical shift in nuclear magnetic resonance (NMR) spectroscopy, or whether the property belongs to a certain class (in particular, the possible existence or non-existence of an object). From the mathematical point of view, we can differentiate several cases of functional relationship between the output training and input data .

Linear regression. The model function is linear and the output value is an experimental value. The output value is obtained as scalar product of the input vector and the weight vector w with the basic equation

Nonlinear regression. The model function is nonlinear and the output value is an experimental value. There is no basic form of an equation, as nonlinear equations can take many different forms. A few of the most flexible curve-fitting forms are n degrees (n > 1), trigonometric functions, exponents and so on:

Classification. The model is nonlinear and the output value is the membership of a class. In this case the non-equality is used (δP is normally introduced to avoid a trivial solution, when all the variables are zero, i.e. the solution is the zero vector):

Depending on the type of function (equality, non-equality, linear or not) several methods can be used to fit the experimental data and parameterize the model. Recursive fitting, regression and classification are most often used. Which technique to use, classification or regression, depends on the results that we want to achieve. If we are trying to predict which class (or category) something belongs to, we use classification. If we are trying to predict a certain value, regression should be used. The difference between classification and regression is defined by the loss function. In regression and classification, the goal is to minimize the value of the loss function. In regression, the loss function increases the further you are away from the correct value for a training instance. In classification, typically, the loss function is 0 for the right class and 1 for the wrong class (perception). The result is that classification only `cares' if it gets the right class, and there is no measure of how `close' it was (or how likely it is that the class is correct). The `goal' of regression is to predict the training values as closely as possible. In regression, different loss functions affect how being close and being really far off are weighted in relation to all training examples. Regression will be also used in Section 6, where structure–property relationships, together with the methods of PCA and dimensionality reduction, will be discussed.

The simplest case of data fitting is `simple linear regression'. Simple linear regression fits a straight line through the set of k points in such a way that it makes the sum of squared residuals of the model (that is, vertical distances between the points of the data set and the fitted line) as small as possible. The solution of the problem is straightforward and well known.

The most widely known application of this tool in crystallography is crystal density prediction. Here, the experimentally determined values of the crystal density serve as the target values and the molecular sum formula serves as the input vector. The sum formula can be considered as the most important input vector in chemistry. The descriptors in the case of a chemical formula are the elements and the input vector contains the frequency of the elements. Many approaches to quantitative structure–property relation (QSPR) and quantitative structure–activity relation (QSAR) analysis use this vector together with some additional descriptors. In our examples we are using it for the prediction of the densities of crystals.

In principle, if a molecule is built up from different sorts of atoms, this is a multiple linear regression task. By assuming that all non-hydrogen atoms of a molecule are of only one sort, Kempster & Lipson (1972) simplified the case to a simple linear regression solution. The result was the so-called 18 ų rule: , where V^sim is the estimated molecular volume, 18 is the weight factor (i.e. w = 18) and is the input vector, which in this case is the sum of all the non-hydrogen atoms in the molecule.

Obviously, this model is very approximate and can be improved by reverting to the `multiple linear regression' task. Multiple linear regression fits a hyperplane through the set of k points in such a way that it makes the sum of squared distances (see Section 3.1) of the model (i.e., the distances from the points to the plane) as small as possible. The equation system that has to be solved is an over-determined linear equation system. The optimal solution of such a system can be obtained with SVD. In this case, the parameterization of a model results in an increment system. The use of increment systems is very simple in practical application. The user just has to count some descriptors, multiply them by the weights and sum them up, which is easy to do. Therefore, increment systems are very popular in crystallography. One of the first increment systems, to predict the sublimation energy of crystals, was developed by Bondi as early as 1963 (Bondi, 1963). Another example is the estimation of the hydration energy of drug-like molecules by counting functional groups in a specific drug set (Pèpe et al., 2002).

The density prediction method of Kempster & Lipson (1972) has been significantly improved by extension of the single descriptor for the non-hydrogen atoms to the sum formula itself (Hofmann, 2002). In this approach, every atom sort i, where i is the atomic number, gets its own number of atoms ν_i and its own `weight' w_i. For simple linear regressionwhereas for multiple linear regression(where 108 is the total number of accessible atom types from more than 110 known at the time of writing) and the density is given bywhere m_molecule is the molecular mass.

In Fig. 14, the results obtained with simple linear regression (red) and with multiple linear regression (green) are compared. It is clear that the predictability, in this case, was improved significantly. The simple but effective approach of multiple linear regression to estimate the molecular volume is widely used these days for cleansing erroneous data from data sets. For instance, the estimate of the density is obtained by multiplying the input vector of the elements by the weight vector w_i, which is the average volume for the element [equation (14)]. A more complex case is that of force fields, which will be discussed in following sections.

Figure 14 Comparing the results of crystal-density predictions obtained with simple linear regression (red) and multiple linear regression (green) for 10 000 crystal structures taken at random from the CSD.

In the common case of linear models the weight vector can be presented as a weight matrix. If the output vector consists of a single value, the weight matrix shrinks to a weight vector. Note that in statistics the application of a `weight' (or `weight vector') is a mathematical device used when performing a sum, integral or average to give to some descriptors more weight, or influence, on the result compared with other descriptors in the same set. Very often, the weight is a simple factor and the score for a possible event is obtained by multiplying the input vector by the weight factor. However, in science the weight factor very often has a physical interpretation.

Nonlinear regression is a form of regression analysis in which observational data are modelled by a function that is a nonlinear combination of the model parameters and depends on one or more independent variables. This kind of analysis is very important in crystal engineering, because along with the classification suggested earlier (Apostolakis et al., 2001), it allows the derivation of intermolecular potentials and their parameterization. Since obtaining the direct solutions of the systems becomes very difficult and time consuming, the data are normally fitted by successive approximations. The `discretization' and `linearization' techniques are used the most often to solve nonlinear problems.

4.3.1. Discretization and linearization

In mathematics, the term `discretization' concerns the process of transferring continuous functions, models and equations into discrete counterparts. This process is usually carried out as the first step towards making them suitable for numerical evaluation and implementation on digital computers. In the simplest case, a continuous function f(r) with (i.e. a real number) is mapped to a discrete function f(k) defined only for the natural numbers , which are the grid points. The set of grid points g(k) can then be chosen as equidistant [for example, all natural numbers ], or obtained with a special function for computational convenience or for different significance in the different regions of the function. For instance, for intermolecular interactions a useful partitioning might be where N_mesh is the mesh size of the interval in which the potential function will be defined. (For example, N_mesh = 200 defines an interval of 10 Å for k = 2.)

The application of this method in knowledge discovery requires balancing the mesh size (the total number of mesh points n_mesh) with the available data. The predictions become more accurate with a larger number of available data, not only in terms of statistical significance, but also because the mesh can be finer.

The term `linearization' in mathematics refers to finding the linear approximation of a function at a given point. In the first step the function has to be discretized and thereafter the function should be linearly interpolated between the mesh points.

To illustrate, the discretization of a Lennard-Jones potential (where r are interatomic distances) iswhere A and B are parameters, and

In Fig. 15, a typical atom pair function for intermolecular interactions is shown. The potential is defined on the mesh points 200/k. At short distances the points are dense, to describe the interaction potential accurately. At large distances they are sparse to economize on the number of descriptors. In between the mesh points the potential is linearly interpolated.

Figure 15 An atom pair potential, which is discretized and linearly interpolated between the mesh points.

4.3.3. The basic functional form of force fields

The basic functional form of a force field encapsulates both bonded terms, relating to atoms that are linked by covalent bonds, and non-bonded `non-covalent' terms. For the non-covalent terms, a large number of functionals exist, which differ in their atom types and the functional dependency of the energy on the distance. One of the earliest potentials used was derived from theoretical considerations about the dispersion and electrostatic interactions. For the long-range interactions it followed an r<sup>−6</sup> law and for the electrostatic term it varied as r<sup>−1</sup>. This potential was then completed with a term for short-range repulsive interactions. For computational convenience an r<sup>−12</sup> (Lennard-Jones, 1931) termor an exponential term (Buckingham, 1938)were chosen. By the 1970s it was recognized that this functional form does not describe hydrogen bonds well, and an additional `10–12 potential' was proposed for hydrogen bonds (McGuire et al., 1972):

Since then, more and more complex forms of potentials have been developed. For example, the shifted force n–m potential (Clarke et al., 1986) iswith The various parameters of the three formulas above can be interpreted by curve sketching or by a machine-learning approach. In terms of curve sketching, ɛ_vdW (E₀) defines the minimum depth of the potential curve (usually referred as a dispersion energy), r₀ is the position of the minimum and σ is the distance at which the value of the potential curve is zero; r_cutoff is introduced for a numerical reason and defines the distance at which the curve should go to zero for all values above the cutoff. A, B and C are parameters that can be interpreted in the terms of machine learning as weight factors for the strengths of the different interactions (see Table 1).

Table 1 The physical significance of the various negative powers Power Interaction Magnitude (kJ mol−1) Name of potential Temperature dependence 1 charge–charge 250 Coulomb   2 charge–dipole 15   T−1 3 dipole–dipole 2   T−1 4 dipole–quadrupole < 1   T−1 5 quadrupole–quadrupole   T−1 6 induced dipole–induced dipole 1 London dispersion   10 hydrogen bond       12 electron–electron   Exchange energy, Pauli energy

4.3.4. Derivation of intermolecular potentials

The method described above in Section 4.3.2 allows derivation of the effective non-covalent interaction without assuming a functional dependency from the outset. The derivation and parameterization of force fields with regression is a more complex case than density estimation. In this case, as a very special representation of the input vector we consider the radial distribution function, RDF. The radial distribution function (or pair correlation function) g(r) in a system of particles (atoms, molecules, colloids and so on) describes how their density (occurrence), ρ, varies as a function of distance from a reference particle. If a given particle is taken to be at the origin O, and if ρ₀ = N/V is the average value of the density of particles, then the density (occurrence) at a distance r from O is . If the descriptors are the discretized frequencies of interactions within the interval r(k) to r(k + 1), then the density ρ₀ is connected with the input vector by

If the RDF is calculated for a single crystal structure, it reflects all the intermolecular distances that occur between all atom types in this crystal structure. The energy of the system, G, can be written as the integral of the RDF with the corresponding atom pair function ɛ. The integral transforms, after discretization, into a sum. For example, in the simplified case of aromatic hydrocarbons with only two kinds of atoms (hydrogen, H, and sp²-hybridized carbon, C) we obtainwhere the weights are the energies of the atom pair function at the mesh points. The weights are still simple factors, but, in contrast to a simple statistical analysis, can be negative as well as positive. The input vectors are obtained from the transformed RDF [equation (25)]. For this case, the number of required potentials results from n_potentials = n_atom types(n_atom types + 1)/2 and for two atom types n_potentials = 3. If the energy is added as a dimension, the problem is a four-dimensional one. The method also works for higher dimensions (we have tested it for up to 67 dimensions), but a four-dimensional space is the limit for reasonable clarity of presentation.

To complete the input vector for deriving the force field for aromatic hydrocarbons, as a first step we extracted from the CSD crystal structures that contain only sp²-hybridized carbon atoms and hydrogen atoms using the `atom sort' filter. This subset of structures was then cleansed. Typical tasks in cleansing are checking for inconsistency in the input vector and checking for its completeness. A further task is the analysis of the data for outliers. For the details and results of this analysis we refer readers to Section 5.

Here the use of filters simplifies the problem of the visualization of the results, as it restricts the number of required descriptors to only two atom types and significantly reduces the dimension of the input vector from roughly 100 to 2. In real applications [see, for example, Hofmann (2010)] we use as a filter the `only organic compounds' option. Obviously, the results obtained using the trained model can only be used for predictions about the corresponding object class: hydrocarbons as in the given case, or organic crystals as in Hofmann (2010).

In total after filtering and cleansing the input vector was composed of 215 structures. On average, these structures had around 15 degrees of freedom and we obtained 3282 equations. The distances for the mesh points were defined at . This choice was advantageous, because it avoided the need to calculate the root of the distances during evaluation of the energy. The experimental data exhibited 90 descriptors; contacts longer than a cutoff of 10 Å were not taken into account and contacts at short distances, below the minimum van der Waals distance, did not occur. The result is shown in Fig. 16. The curves obtained are similar to common interaction potentials, but they are not completely smooth. Since we did not presuppose a functional dependency between the distance and the potential, the potentials show some noise. As a reference, in Fig. 16 we have added the radial distribution functions for the corresponding contacts. To a first approximation one would expect that the minimum of each potential should coincide with the maximum of the RDF and the curves should be related to the Boltzmann law, = , where k_B is the Boltzmann constant. In fact, the minima are shifted to slightly higher distances than the RDFs suggest.

Figure 16 The potentials for aromatic hydrocarbons obtained by singular value decomposition (SVD) are shown in dark red, dark green and dark blue. For comparison, the corresponding radial distribution functions (pale red, pale green and pale blue) are also shown.

It is important to mention that the solution of the equation system using SVD indicated that the matrix was not full rank. It gave a solution for only 88 weights, and two weights could not be determined. This kind of result is common. In general, some equations can appear as co-linear. An obvious reason for this is if several descriptors occur only in a single experimental data point. In this case, the weights of these descriptors can assume arbitrary values and a defined value only exists for the sum of them. For example, this happened for the structure with the CSD refcode SALFUG (Fig. 17), which has an intermolecular H⋯H interaction of length 2.098 Å. This singular distance contributed to the frequency vectors of the two adjacent descriptors: as the floor is the same as the ceiling , these results were co-linear. Only a visual check of the structure allowed determination of whether the structure was correct or whether some inaccuracies in the H-atom position were present, hence this structure had not been removed by cleansing earlier. The molecule is quite unusual, which might be responsible for unusual contacts as well. This anomaly could only be detected at the stage of SVD.

Figure 17 The structure with the CSD refcode SALFUG has a very short H⋯H distance and was recognized as unusual during multiple regression.

Another method for the parameterizing of force-field models is classification. This method, using a support vector machine (SVM) algorithm, was described in detail and used in Hofmann & Apostolakis (2003), Apostolakis et al. (2001), and Hofmann (2010).

In classification, if data can be divided into different classes, then the parameters of a model can be optimized to separate these classes as effectively as possible. Often only two classes are considered, one `true' and one `false'. The advantages of the method are that the parameters can be determined for an arbitrary model, the model can be multidimensional and any arbitrary function for the functional dependency is allowed. The largest drawback of the approach becomes apparent as soon as we begin to consider optimization. Any algorithm for optimization results in a local minimum, and no one algorithm can be guaranteed to find a global minimum. The more complex and multidimensional the input data are, the more likely it is that the algorithm does not find the global minimum. As a result, iterations of the parameterization with different starting points will finish with different values, whereas (and this should be emphasized) parameterization of the force field by regression give us the unique solution. Therefore, force-field parameterization by classification requires a very good initial guess. Nevertheless, the method is very attractive because of its tolerance of any model and its broad applicability. The problem has only two requirements for its application:

(1) for any experimental data a reference state is known, which is higher in its score (energy); and

(2) the experimental data are a local minimum of the model.

If these two requirements are fulfilled, we can always define two classes: `true' and `false'. The first class, which we will call `experimental data' or `true' henceforth, can be defined with the help of the first requirement: it contains the experimental data minus the reference state. The scores of all elements in this class must be zero. In crystals, this means explicitly that the interaction energy within a crystal must be below the energy of the free molecules. If the conformation energy remains constant or there are only trivial changes in the conformation during crystallization, it is equivalent to the statement that the free energy of the crystal formation must be below zero. The second class, which we will call `virtual data' or `false' henceforth, can be defined with the help of the second requirement: it contains the virtual data minus the reference state. The virtual data can be generated in the case of crystal structures by applying some small distortions to experimental structures. From the second requirement it follows now that the score (the free energy) of the distorted structure minus the score of the experimental structure must be above zero. In total, we obtain the class `true', where each element must have a score below zero, and the class `false', in which each element must have a score above zero. In Fig. 18 we show the results of classification for the case of aromatic hydrocarbons, described above in Section 4.3.4. The set of the data to classify was composed of the 215 experimental structures used above plus the set of the virtual structures (around 10 virtual structures are generated for every experimental structure). The algorithms that allow the parameterization to achieve an optimal separation of the two classes are SVMs. They were first used in crystallography as early as 2001 (Apostolakis et al., 2001).

Figure 18 Results of the classification according to the energy (vertical axis, interval 200 score points). The structures on the horizontal axis are ordered according to their scores (from the lowest to the highest value). All the experimental structures (blue) have energies below zero and most of virtual structures have energies above zero (green). However, some virtual structures have energy values below zero (red).

As pointed out above, the first step of the method is to define a good starting guess. As a starting point, the parameters of one of the common force fields or the force field obtained by regression (Section 4.3.4) can be used. The possibility of using the force field obtained from regression opens new horizons in developing highly customized potential functions for special groups of compounds, or for compounds under unusual or extreme conditions. A general possibility for overcoming the problem of the unknown functional dependency of the atom pair interaction as function of the distance is to develop the function as a series. Obviously, the terms of the series are forbidden to be co-linear and they should fulfil the limit conditions. In the case of atom pair functions, one implication of this is that for infinite distance the interaction energy has to be zero. It is possible to use reciprocal power series of the form , with n ≥ 2. With this condition for n the common problem of the divergence of the Coulomb term is avoided and the Coulomb energy during parameterization becomes intrinsically included into the other terms. If n = 6 is chosen and the last term is taken as r⁻¹², the potential corresponds to the common Lennard-Jones potential. A full list of the different terms and their physical meaning can be found in Table 1.

Here we demonstrate the method for the aromatic hydrocarbons, the same example as used in Section 4.3.4. In the first step we generated the false (virtual) structures, which we then used for training. The number of searched atom pair potentials was three, as before. Series of the form were used. Any such series contains four terms and, in total, 12 parameters have to be determined. The initial guesses for the parameters were obtained by fitting the functions to the atom pair functions obtained by regression in Section 4.3.4. Starting with this fit, the sum of the distances of the outliers from the hyperplane was minimized by a simplex minimization. In such a procedure, the distance for the experimental structure is given by , while for the virtual data the error is given by . Since the system has a trivial solution by setting all parameters to zero, a boundary condition needed to be introduced during minimization. The functions were always scaled to keep the minimum value at −10; .

4.4.1. Outliers in classification (energy cleanser)

Special attention should be paid to the analysis of the results of classification. In Fig. 19, we plot the energies of the virtual structures against their occurrence in the interval from −0.1 to 0.1 score points, a total range of 0.2 score points, enlarged from Fig. 18.

Figure 19 More detailed results of the classification according to the energy (vertical axis). The structures on the horizontal axis are ordered according to their scores. Only the virtual structures are shown. The energies of most of the virtual structures are above zero (green). However, some virtual structures (outliers) have values below zero (red).

As one could see in Fig. 18, all the experimental structures have an energy below zero and most of the virtual structures have energies above zero. However, and this can be seen on the enlargement shown in Fig. 19, some virtual structures have energy values below zero (red points). The fact that all the experimental structures are assigned to negative energy shows that the internal cleansers of the CSD work very well and, at least for the present set of training structures, no structure with catastrophic errors was found. This was not always the case (Hofmann & Lengauer, 1997), and we should thank the continuous improvements in the data cleansing during the last 20 years (Groom et al., 2016) for the good quality of the present data.

The importance of outlier analysis for model development and data cleansing was discussed in Hofmann & Kuleshova (2010). Here we give some further insights to the problem. As we have mentioned above, most of the outliers were detected in the virtual data. It was noted that 1.89% of the virtual structures have energies below zero. This means that some distorted structures have a lower energy than the experimental structures.

This can happen if the model is not accurate enough and/or (for example) too few parameters (descriptors) are available. Another possible reason can be errors in the experimental data. Therefore, it is always worth looking at the outliers in detail to work out whether it is possible to improve the method. In this case, the strongest outlier was the structure ANNULE10 (the isolated red point on Fig. 19). This single structure contributes 12% of the overall error. Fortunately, the database contains a second entry for the same compound, with the CSD refcode ANNULE01. The direct overlaying of the molecules, Fig. 20, from the two structures immediately reveals a problem: in ANNULE10 (green), some of the hydrogen atoms are positioned inaccurately. The C—H bond is a little bit too long and the hydrogen-atom position is too far out of the plane of the sp²-hybridized carbon atoms. The length of the C—H bond can be normalized easily by using standard bond lengths for hydrogen in the programs Mercury (Macrae et al., 2020), FlexCryst (Hofmann et al., 2016) and Materials Studio (Accelrys Software Inc., 2013). The position of an `out-of-plane' hydrogen atom cannot be corrected automatically yet. As another result of the bad positioning of the hydrogen atoms, one H⋯C interaction is shortened from 2.803 Å to 2.599 Å. These small errors in hydrogen-atom positions were detected only during the classification, and cause the structure ANNULE10 to be the largest outlier.

Figure 20 The overlay of molecular fragments in ANNULE01 (grey) and ANNULE10 (green) illustrates the error found by the classification. ANNULE10 has a C⋯H contact that is too short: 2.599 Å.

To our knowledge, systematic errors are not common in databases, since the error can usually be eliminated if the cause of the systematic error can be identified. The methods for detecting a systematic error are well known in statistics. A simple check is whether the measured average values coincide with the expected average values or not.

Conceptually, the detection of systematic errors during supervised classification is not completely automatic; nevertheless, it might be the most powerful method. What commonly happens in the case of systematic errors is that during supervised classification a large part of the data cannot be assigned to the correct class and stay as outliers. In the next step these outliers have to be analysed and the crystallographer has to use their intuition to find a descriptor that differentiates between the correctly classified data and the outliers. If such a descriptor is identified, the origin of the error can be traced back.

Density cleanser. The effectiveness of this approach can be illustrated by the inclusion of the crystal density as an additional descriptor for monitoring crystal structures. For example, during the classification of the crystal structures with space group P1 performed by Hofmann (2002) (according to the scheme described in Section 4.4) a huge number of outliers were observed. Further analysis of the two classes obtained with linear regression for the crystal density (Section 4.2) revealed different weights for the average volume of the atoms in these groups. The values for the major class were less than half of those for the minor class (the outliers). Attempts to trace a possible reason for this suggested that in some cases an overbar had been lost during text processing, and that the structures in the minor class in fact belonged to the space group . In the past, several systematic errors have been detected; other examples are errors in the setting (for example P2₁/a ↔ P2₁/n ↔ P2₁/c) or the handedness in certain space groups (P2₁2₁2₁). We have already mentioned that a density cleanser now is one of the cleansers used by the CSD. The other example of anomaly detection with classification was described in Section 4.4.1 and was referred to as an `energy cleanser'.

Random errors are always present in experimental data. They depend on the (possible) accuracy of the data acquisition. They are easily recognizable, within the methods of data mining, because the data are scattered around the line in regression, around the hyperplane in multiple regression, or the clusters are blurred during classification. But the errors do not show any structure. Only the estimation of the average value can be improved by cleansing (preprocessing) the experimental data set.

A single set of data can be only improved by changing the data acquisition. A common problem in this context is the determination of hydrogen-atom positions by X-ray crystallography. Since this method detects the electron density rather than the positions of the nuclei, the positions of the nuclei of atoms without electrons cannot be determined. Therefore, the positions of highly acidic hydrogen atoms are very often ill-conditioned. The only ways to overcome the problem are to change the method and collect data using neutron diffraction, or to set the relevant values to average ones (i.e. use normalized hydrogen-atom positions).

RDF cleanser: That some atoms are ill-positioned becomes obvious if the radial distribution functions (RDFs) are plotted as in Fig. 21. It is seen that at short distances (and especially for hydrogen atoms) some contacts fall outside normality and form isolated clusters.

Figure 21 Average RDFs for the set of aromatic hydrocarbons (shown in dark red, dark blue and dark green) used in Section 4.3.4, and the RDFs for a single structure (pale red, pale blue and pale red), CSD refcode BOXGAW.

For the detailed investigation of the outliers in this case, a measurement of anomaly (similarity η) has to be introduced. This index has to be high for the outliers and low for normal structures. For this aim the relative difference between the RDF of the sum of the structures, RDF^set, and the RDF of a specific structure, RDF^structure, can be used [equation (27)]. A good measure for the value of η for the anomaly will be the maximum of this relation at the mesh points. Furthermore, it is desirable that the error should have a Gaussian distribution, since Gaussian distributions can be analysed with descriptive statistics, which are highly developed. For instance, they allow us to differentiate between systematic and random errors and to estimate the likelihood that the structure is wrong. Ratios can be mapped to Gaussian distributions by taking the logarithm, which gives us an equation for differentiating between ill- and correctly positioned atoms:

The result can be plotted as histogram. In Fig. 22 the result for all organic crystal structures containing common atom types is shown. The set consists of 74 625 crystal structures. The anomaly of the structures rises with increasing value of η, while with decreasing η the structures become more and more similar to the average. At the mean the structures have a score of around two (mean value of η: 2.05; median: 1.97).

Figure 22 Anomalous structures obtained with the RDF cleanser.

The most anomalous was found to be the crystal structure VEXXUS (highlighted in the histogram by the red arrow). VEXXUS (Gallardo et al., 2007) is a cocrystal of tetraethylammonium with 1,3,5-trinitrobenzene and, in fact, possesses a very particular crystal structure: analysis revealed very unusual contacts between the molecules of VEXXUS (see Fig. 23), including a interaction with a very short distance of 2.478 Å.

Figure 23 According to the RDF cleanser results shown in Fig. 22, the most anomalous structure is VEXXUS, with a very short C⋯C distance of 2.478 Å.

The ten most anomalous structures are listed in Table 2. They have errors above 9.85.

Table 2 The most anomalous structures according to the RDF cleanser CSD refcode Similarity η Interaction type r (Å) VEXXUS 11.2 2.50 WABKAL 10.7 2.62 IVEZAK 10.3 2.62 PARNIG 10.2 C—H⋯ 1.96 VUCFUU 10.2 2.32 UXIPEX 10.1 2.54 SEQBUN 10.0 C—H⋯ 2.02 PYRIDO 10.0 2.62 LAJWIE 9.9 2.88 VAMVEI 9.9 C—H 1.94

Model errors can only be detected during supervised learning. They can be recognized when a lot of data cannot be assigned to the right class or have bad correlation coefficients during regression and the output data can be reflected only poorly by the trained model. To understand that the data are reflected poorly, it is necessary to examine the output data during machine learning, and in this sense only supervised learning is able to detect whether some important descriptors are missing. An obvious strategy for solving this problem is to extend the number of descriptors and to collect all data for any descriptors that can be thought of. Below we highlight the importance of the proper choice of descriptors when looking for errors in models and missing descriptors.

In Section 3.3 on the clustering of crystal structures, it was demonstrated that clustering allows us to recognize descriptors that have been overlooked, like the temperature and pressure. It is astonishing in this context that methods of CSP take no account of the effects of temperature and pressure and yet are successful (Reilly et al., 2016), even though some experimental crystallographers specialize in high-pressure or low-temperature studies. For example, crystallography under high pressure has significance for the deep sea (for instance, in the handling of the oil spill in the Gulf of Mexico), whereas crystallography at low temperature is relevant to space and astronomy.

In Section 4.2 it was shown that prediction of the crystal density works much better within the multiple linear regression when a weight factor is added to the single descriptor of the number of atoms.

In Section 4.4, however, when discussing the analysis of intermolecular interactions in crystals it was demonstrated that extending the descriptors has limits and that the number of descriptors cannot become infinite. The number of descriptors must always be lower than the number of data, otherwise the problem gives rise to an undetermined equation system. If the connectivity between the input data and the output data is a function, this limits the mesh size that can be applied when discretizing the function. The accuracy increases when a finer mesh is chosen, but the mesh size is limited by the available data.

In Section 6.1, we will show how to reduce an infinite number of descriptors to a reasonable number of descriptors by principal component analysis, PCA, for the prediction of crystal solubility.

Currently the use of PCA is the poor cousin in machine learning. Still, computing power is growing quickly enough to keep pace with the amount of data being collected. While small centres, like the Center for Advanced Studies, Research and Development in Sardinia (CRS4), can store some petabytes, the largest data storage centre of the world, the Range International Information Group data centre in Langfang, China, stores some thousands of petabytes. A theoretical limit for the present technology is that common processors with a 64-bit address bus can address a maximum data set with 16 000 petabytes of entries.

Principal component analysis (PCA) is probably the most popular multivariate statistical technique and it is used by almost all scientific disciplines. It is also likely to be the oldest multivariate technique. It can be traced back to Pearson (1901) and for more details we direct the reader to, for example, Abdi & Williams (2010) and Samudrala et al. (2014). PCA analyses a data table representing observations described by several dependent variables, which are, in general, intercorrelated. Its goal is to extract the important information from the data table and to express this information as a set of new orthogonal variables called principal components. PCA also represents the pattern of similarity of the observations and the variables by displaying them as points in maps (Jolliffe, 2002).

The main linear technique for dimensionality reduction, PCA, performs a linear mapping of the data to a lower-dimensional space in such a way that the variance of the data in the low-dimensional representation is maximized. In practice, the correlation matrix of the data is constructed and the eigenvectors of this matrix are computed. The eigenvectors that correspond to the largest eigenvalues (the principal components) can now be used to reconstruct a large fraction of the variance of the original data. Moreover, the first few eigenvectors can often be interpreted in terms of the large-scale physical behaviour of the system. The original space (with dimension equal to the number of points) is reduced (with data loss, but hopefully retaining the most important variance) to the space spanned by a few eigenvectors.

6.1.1. Molecular superposition

PCA is usually explained via an eigen decomposition of the covariance matrix. However, it can also be performed via singular value decomposition, SVD, of the data matrix. An illustrative example of the relationship between SVD and PCA, with which many readers might be familiar, is the application of SVD to the positions of the atoms in a molecule. In Euclidean space the atoms have three descriptors for their position: x, y, z. With SVD we obtain three weight vectors. The first gives the direction in which the molecule is most extended; the second gives the width and the last gives the flatness. The singular values give the variances for the three directions. The more dispersed the points are, the larger the singular values. If the molecule is completely flat one singular value is zero and one dimension is redundant. If the molecule is linear, two singular values are zero and the dimensionality can be reduced to a single descriptor. Let us apply SVD to the determination of the axis of inertia in a molecule. The matrix of interest in this context is the tensor of inertia. In contrast to the more general matrix in equation (18), which is rectangular, in the case of the tensor of inertia the matrix shrinks to a quadratic 3 × 3 matrix. In this case SVD becomes equal to solving an eigenvalue problem. The weight vectors are identical to the eigenvectors and the eigenvalues correspond to the singular values. Applying SVD gives as a result the three weight vectors corresponding to the axes of inertia. The largest of the singular values is the moment of inertia for the principal axis. Using this approach gives us a very quick and efficient algorithm for the superposition of molecules. For two molecules to be compared the principal axis of each has to be determined and a principal-axis transformation performed. Then a similarity index can be calculated, for instance, the root-mean square of the shortest distance between the atoms in the molecules. An advantage of this algorithm in comparison with other algorithms is that the molecules do not need to have the same number of atoms. An example, the superposition of biphenyl and fluorene molecules, is shown in Fig. 24.

Figure 24 Superimposing biphenyl and fluorene with singular value decomposition. The calculated inertia ellipsoids are a measurement for the similarity. For flat molecules one dimension becomes redundant.

In the previous sections many different sub-areas of machine learning have been described, and the reader might get the impression that everything can be easily solved by machine learning. In this section another method will be presented. One of the main aims of physicists and other scientists is to find out experimentally how an event (the output data) depends on different influences and to try to exclude unimportant influences. In this sense the experimental work is an experimental PCA. Therefore, we can simplify a problem a great deal by taking well known results into account. The example chosen here is the pharmaceutical application of crystallography. Most drugs are microcrystals, and crystallography plays an important role in the search for new drugs, their patenting and in quality control. With the arrival of combinatorial chemistry, it has become possible to prepare millions of compounds in a single process. This raises the question: Is it possible to predict the solubilities or biological activities of the compounds, or it is simpler to determine them by high-throughput screening?

The natural approach using data mining is to describe the molecules as much as possible using descriptors, perform PCA, eliminate the redundant descriptors, and to train the model by a supervised method, as described in Section 4. This scheme can be used to predict an arbitrary property of a molecule or a crystal. Schemes that use multiple linear regression (Sections 4.2 and 4.3.4) are the most popular. In this case, the result is an increment system, which is easy to use. If the aim is to predict an arbitrary property the scheme is called a quantitative structure–property relationship (QSPR) method, and in the specific case of biological activity the scheme is called a quantitative structure–activity relationship (QSAR) method. An early example of QSPR in chemistry was the prediction of the melting points of alkanes in 1842 (Kopp, 1842). Recent applications in inorganic crystallography can be found, for example, in Rajan (2010), Samudrala et al. (2014), Villar et al. (2008) and Rajan et al. (2009).

The accuracy of this method grows with the increasing numbers of descriptors. In the early days, experimental data were used as input data. Since the measuring of input data for large data sets (such as from combinatorial chemistry) can become tedious, these days calculated descriptors are very often used. To enable the use of as many descriptors as possible, programs have been equipped with estimates of more and more descriptors (Todeschini & Consonni, 2009). One of the leading programs in this field, Dragon, calculates 5270 molecular descriptors at present (Mauri et al., 2006). The list of descriptors includes the simplest (atom types, functional groups and fragment counts), as well as topological and geometrical descriptors, three-dimensional descriptors, quantum-mechanical descriptors, properties and so on.

One of the first to perform an intuitive reduction of this flood of descriptors for QSAR methods was Lipinski (Lipinski et al., 2012). He ended up with `the rule of five'. Only four descriptors were assumed to be valuable. In addition, the input vectors are logicals (true or false), and are weighted equal (true). The simplicity of equation (28) makes this rule one of the most popular (the article has been cited more than 10 000 times).

Looking at the problem of biological activity of oral drugs in more detail, it becomes clear that it can be divided into several aspects. The oral drugs must go into solution and then penetrate the cell membranes. Consequently, methods exist that try to predict separately the solubility or the intestinal permeability. To predict the solubility of an innovative drug at a very early stage, ideally before synthesis, is a challenging problem of high pharmaceutical significance. Fortunately, the solubility is easy and accurate to measure. Therefore, large and accurate databases are available, for instance at the National Institute of Standards and Technology (NIST). This makes the idea of applying QSPR to this question very attractive, and a large number of studies of this kind have been undertaken.

The large number of descriptors in QSAR methods is in contrast to the presentation given in many textbooks. In thermodynamics, the solubility, K_s, of a crystalline substance, for the case of simple dissolution, depends mainly on two kinds of interactions. The first defines how tightly the molecule is bound to its own crystal lattice, G_cryst, and the second how strongly the molecule associates with the solvent, G_solution (or G_hydration for the solubility in water). In thermodynamic equilibrium these two energies have to be equal. The solubility itself is a function of the solution energy in an infinitely dilute solution (the `standard energy of solution') and the concentration. The free energy of solution depends in an additive way on the standard energy of solution and the concentration:

The free energy for a solution at a particular concentration is connected, to a first approximation, to the concentration by Boltzmann's entropy formula, which connects the entropy for a system with the energy. The considerations are based on the partition function for a number of equal particles, therefore the concentration has to be expressed in terms of mole fractions, X. Replacing the free energy for a solution at a particular concentration we obtainwhere R is the universal gas constant and T is the temperature. can be expressed as a linear regression: where a and b are parameters.

It is problematic, however, that neither G_cryst nor G_solution is easily accessible by experiment or by simulations. Nevertheless, recent progress in crystal structure prediction (CSP) allows us to estimate the free lattice energy of the crystals, G_cryst, with passable computational effort. Even if crystal structure prediction by itself still requires enormous computational effort, the lattice energy can be quickly estimated by a simplified version of CSP. It is important in this context that during CSP a lot of polymorphs are found within a very narrow range of energy. For predicting the solubility, it is sufficient to estimate the lattice energy roughly and knowledge of the real crystal structure is not needed. In Fig. 25 the calculated lattice energies of some experimentally determined crystal structures and the lattice energies of their top-ranked predicted structures are compared. With this simplified CSP the lattice energy is constantly predicted as a little bit lower than the real lattice energy, but the correlation is nearly perfect, as the correlation coefficient of 0.99 shows. (A list of the drugs and drug-like compounds that feature in Figs. 25–27 is given in Table 3.)

Figure 25 Prediction of the lattice energy (in kJ mol−1) with FlexCryst.

The other descriptor, i.e., the solution or the hydration energy, can be obtained from an increment system (Pèpe et al., 2002) or by molecular dynamics. It is quicker to evaluate the hydration energy using the increment system; however, using molecular dynamics is more precise. The increment system is based on 2D descriptors of the molecule, i.e. the input vector contains the frequency of different functional groups of the molecule. Combining both descriptors and performing a two-dimensional regression according to equation (31) gives a good correlation between the experimental and the predicted solubility.

For further improvement of the prediction more descriptors can be added to the formula. One of the most common descriptors is the molecular volume (Lüder et al., 2009). The result of three-dimensional regression with this additional descriptor is plotted in Fig. 26. The correlation between the solubility and the molecular volume is negligible (the correlation coefficient of 0.01 is tiny). This means that when only two descriptors – the lattice energy and the hydration energy – are used the molecular volume becomes redundant.

Figure 26 The correlation between solubility and molecular volume. The molecular volume is redundant.

A comparison of the predicted solubility method with other methods is shown in Table 4. The data are taken from Duchowicz et al. (2008) and show the success of the different methods for a set of 21 drugs. The r.m.s.d. value of 0.85 is slightly higher than for the best model of Hou et al. (2004) (r.m.s.d. = 0.664) and the model of Huuskonen (2000) (r.m.s.d. = 0.810). However, their models both use more descriptors than the size of the training set and the ratio of the number of parameters to the number of data sets is 4:1 for the model of Hou et al. (2004) and 1.5:1 for the model of Huuskonen (2000). In both cases this is questionable; overfitting might have occurred. The method presented here is better as it uses only two descriptors, and the ratio of the number of parameters to the number of data sets is much lower at 1:17.

In Fig. 27 the correlation between the predicted and experimental solubility is shown. The experimental data were taken from the online database of the Virtual Computational Chemistry Laboratory (2016). The data set covers a broad range of solubility and varies over seven orders of magnitude. Fenbufen, for example, is nearly insoluble, with a pK_s of −5.06, in contrast to acetamide, which is nearly mixable with water (pK_s = 1.58 ≃ 2.25 kg l⁻¹). The correlation coefficient of 0.80 compares well with other methods. By PCA it is possible to show that the irrelevant descriptor, the molecular volume, cannot improve the prediction. As suggested in many textbooks, it is very likely that the concentration should enter in some way as a descriptor. It is well known that for high concentrations, as for urea and acetamide, the hydration shell around the molecules is incomplete because of a lack of water, and that the solute molecules interact directly with each other. To correct this problem activity coefficients are added to the first approximation of equation (30).

Figure 27 Solubility prediction with reduced dimensionality (only hydration and lattice energies considered).

Another example, and one of how to estimate the solubility even if some very important data are missing, is the case of cocrystals (Kuleshova et al., 2013). In pharmacy, cocrystal formulation is very common, with the aim of improving the solubility of the active pharmaceutical ingredient (API), so it would be very valuable to be able to estimate, even if not the solubility by itself, then the change in solubility after cocrystal (CC) formation with one of the possible coformers (CFs). For crystals of an API, the basic expression for solubility can be written asHere, as in Section 6.2, the first term, , defines how tightly the molecule is bound to its own crystal lattice and the second, , how strongly the molecule associates with the solvent. Within the approach described above (see Sections 4.3.4 and 4.4), the calculation of the crystal energy of any molecular system is not a problem, and takes very little computational power.

In recent years several other relatively simple and quick approaches have been developed (Karamertzanis et al., 2009; Musumeci et al., 2011; Perlovich & Raevsky, 2010). On the other hand, molecular dynamics (MD) estimation of hydration energy requires complex simulations and high computational power (Orozco & Luque, 2000; Wan et al., 2004; Nagy, 2004; Duffy & Jorgensen, 2000; Åberg et al., 2004; Jorgensen & Duffy, 2000).

A simple mathematical reformulation allows us to exclude the second term and to avoid the tedious calculations of a hydration energy.

Let us assume that an API can form different cocrystals (CCs) with a variety of coformers (CFs). The simple dissolution of a cocrystal in water, Fig. 28, is described by the reaction with the equilibrium constant K_s.

Figure 28 Scheme of simple dissolution of a cocrystal in water.

For the case of a two-component 1:1 cocrystal, equation (32) will transform to

As soon as we need just information on the change in the solubility, it is sufficient to consider the difference between equations (32) and (34), which will characterize the `relative solubility' and, in addition, will allow us to exclude from consideration the free energy of hydration of the API, :

Further simplification of the equation is achieved by replacing the free energy of hydration of the CF, , by After consequent transformations the final expression for the relative solubility, , can be written as where From expression (37) it follows that the relative solubility of the cocrystal will be defined by the solubility of the coformer, CF, and by the value of ΔG calculated according to equation (38). For cases where the solubility of the API and coformers are known, the solubility of the cocrystals can be calculated directly:

Data for pharmaceutically acceptable coformers are often stored in specific libraries and their physico-chemical characteristics are known. In these cases expression (37) allows a quick estimate of the relative solubility (virtual screening) of extensive numbers of practicable cocrystals.

To illustrate how this scheme works, we found the structures of four cocrystals (COLHER, WIKTEP, DAQYOK and RILJEB) and their pure components with known data on solubility in the CSD (Table 5). The APIs are the flavonoid naringenin (DOLRIF) and 4-aminobenzoic acid (AMBNAC); as coformers, CFs, we have triethyldiamine (TETDAM), trinitrobenzene (TNBENZ), 4-nitrobenzoic acid (NBZOAC) and 4-nitroaniline (NANILI). For all of them crystal structure predictions were performed in accordance with the protocol described in Section 3.4 and the values of ΔG were calculated according to equation (38).

Table 5 The estimation of the relative solubility for some cocrystals Cocrystal (mol l−1) (mol l−1) ΔG (kJ mol−1) (mol l−1) COLHER (DOLRIF:TETDAM) 0.002 4.012 4.690 22.5 0.045 WIKTEP (AMBNAC:TNBENZ) 0.398 0.001 6.758 0.153 0.061 DAQYOK (AMBNAC:NBZOAC) 0.398 0.002 3.500 0.007 0.003 RILJEB (AMBNAC:NANILI) 0.398 0.004 2.530 0.010 0.004

In all cases ΔG was found to be positive, which means that our calculations suggest that cocrystal formation is favourable. As one can see, the cocrystallization of naringenin with much more soluble TETDAM helps to increase the solubility of the flavonoid by more than 20-fold. Crystallization of AMBNAC with the poorly soluble coformers reduces the solubility of the API. In these cases the values of the relative solubility are less than unity. For coformers with comparable solubility (in our example TNBENZ, NBZOAC and NANILI) the values of ΔG begin to play a leading role in the final solubility of the cocrystal, so the solubility of WIKTEP is higher, owing to the higher value of ΔG.

Further developments will obviously be made possible by growth in both databases and computer power. These will allow for the introduction of additional descriptors, which are still considered to be beyond the current scope of data mining. For instance, data-mining force fields only recently reached an accuracy that allows them to consider the temperature and, accordingly, to describe thermal effects in molecular crystals (Hofmann & Kuleshova, 2018). This has become accessible since data for variable-temperature experiments were added as a separate descriptor to the CIF template used by the CSD. The inclusion of more descriptors in the CIF template (for example, pressure, sublimation energy and so on) would open new possibilities for knowledge discovery.

Data miners are also currently developing methods for combining data from various databases. Combining existing databases is problematic for commercial reasons, but more and more algorithms are appearing that allow data to be used from a combination of various databases. This is not only the case for databases of the same class, such as the structural databases (the CSD, PDB and ICSD); it is also possible to combine data from databases of different classes, such as from thermodynamic databases (for instance, from the NIST database and from structural databases). The latter would allow the scaling of score functions obtained with data mining to standard energy units.

As a growing problem, one should mention the appearance of `fake data' in databases. Because it is very often simpler and cheaper to use an existing model to calculate a specific property rather than to measure it or conduct a literature search for experimental data, people tend to confuse virtual data and experimental data. There is thus a danger that the line between the real and the virtual may be blurred. The organizations behind the databases are well aware of this problem; however, we sometimes find structures in the CSD that have been calculated but not measured.

There is also progress at present in `hybrid methods', where the crystal structure is calculated within a common quantum-chemistry approach and only the exchange correlation term of density functional theory is obtained by data mining.

The further development of algorithms is difficult to predict, but we have presented in this chapter instances in which the weight vectors and the vector matrices have a physical meaning. This should help the reader to understand the methods. However, neural networks, for instance, can be used as a completely `black box' tool. While the trained network might predict perfectly, and better than any other model, this has some drawbacks that contradict the classical approach of science. Because black boxes do not produce any knowledge, they cannot extrapolate, and even interpolation can fail. But one can hope for significant progress with such approaches in the future.