Mapping and Classifying Molecules from
a High-Throughput Structural Database
High-throughput computational materials design promises to greatly accelerate the process of discovering new materials and compounds, and of optimizing their properties. The large databases of structures and properties that result from computational searches, as well as the agglomeration of data of heterogeneous provenance leads to considerable challenges when it comes to navigating the database, representing its structure at a glance, understanding structure-property relations, eliminating duplicates and identifying inconsistencies. Here we present a case study, based on a data set of conformers of amino acids and dipeptides, of how machine-learning techniques can help addressing these issues. We will exploit a recently-developed strategy to define a metric between structures, and use it as the basis of both clustering and dimensionality reduction techniques – showing how these can help reveal structure-property relations, identify outliers and inconsistent structures, and rationalise how perturbations (e.g. binding of ions to the molecule) affect the stability of different conformers.
pacs:Valid PACS appear here
Computational materials design promises to greatly accelerate the discovery of materials and molecules with novel, optimized or custom-tailored properties. With this goal in mind, several community efforts have emerged over the past few years aiida (); clean_en (); esp (); oqmd (); pauling (); matgenome1 (); matgenome2 (); PhysRevLett.108.058301 () that aim at generating, and/or storing large amounts of simulation data in publicly available databasesPhysRevLett.114.105503 (); PhysRevB.92.014106 (); PhysRevB.92.094306 (); kusne15screport (); ramkrsinan_2014sd (); PhysRevB.90.155136 (); PhysRevB.90.155136 (); db-paper (). The development of these repositories of structural data, and of associated materials properties (e.g. formation energy, band gap, polarizability, …) poses considerable challenges, from the points of view of guaranteeing consistency, accuracy and reliability of the stored information, as well as that of extracting intuitive insight onto the behavior of a given class of materials and of data-mining in search of compounds that exhibit the desired properties or that are somehow interesting or unexpected.
In order to automate these tasks – which is necessary to unlock the full potential of computational materials databases that can easily contain millions of distinct structures – a number of different machine-learning algorithms have been developed, or adapted to the specific requirements of this field rodr-laio14science (); clustering-rui (); gang_clustering (); cartography (); prasanna-sd (); ferg+10pnas (); ceri+11pnas (); trib+12pnas (); ceri+13jctc (); rohr+arpc13 (). A fundamental ingredient in all of these approaches is a concise mathematical representation of a molecular or crystalline structure, that can take the form of fingerprints (low-dimensional representation of the structure of the atoms) or more abstract measures of the (dis-)similarity between elements in the database, such as distance or kernel functions.
In the present manuscript we will present a demonstration of how a combination of a very general approach to quantify structural dissimilarity de+16pccp () can be combined with non-linear dimensionality reduction and clustering techniques to address the challenges of navigating a database of molecular conformers, checking its internal consistency and rationalising structure-property relations. Even though we will focus in particular on a energy/structure data set of amino acid and dipeptide conformers obtained by an ab initio structure search aminodb (); db-paper (), many of the observations we will infer are general, and provide insight on the application of machine-learning techniques to the analysis of molecular and materials databases generated by high-throughput computations.
I A Toolbox for Database Analysis
In this section we will describe a specific combination of atomic structure representations and of analysis techniques that can be used to deal effectively with large databases of materials and molecules. We will however also briefly summarize some of the alternative approaches that could be used to substitute different components of our tool chain.
i.1 Fingerprints and Structural Dissimilarity
The most crucial and basic element in any structural analysis algorithm is to introduce a metric to measure (dis)similarity between two atomic configurations. Many options are available, with different levels of complexity and generality, starting from the commonly used Root Mean Square (RMS) distance. In order to deal with symmetry operations or condensed phase structures, several “fingerprint” frameworks have been developed sprint (); PhysRevLett.108.058301 (); PhysRevB.90.104108 (); PhysRevB.89.235411 (); pilania13screport (); PhysRevB.88.054104 (); rupp+07jcim (); hirn+15arxiv (); qm7b (); PhysRevLett.108.253002 (); PhysRevB.92.045131 (); anatoleIJQC (); bag.of.bonds (); newstefan (), that assign a unique vector of order parameters to each molecular or crystalline configuration: a metric can then be easily built by taking some norm of the difference between fingerprint vectors. Any of these distances could be taken as the basis of the classification and mapping algorithms that we will describe in what follows.
In this paper we will use instead a very flexible framework (REMatch-SOAP) that is based on the definition of an environment similarity matrix , which contains the complete information on the pair-wise similarity of the environment of each of the atoms within the molecules and . In our framework, the similarity between two local environments and is computed using the SOAP kernel PhysRevB.90.104108 ()
The REMatch kernel is then defined as the following weighted combination of the elements of
where the optimal combination is obtained by searching over the doubly stochastic matrices the one that minimizes the discrepancy between matching pairs of environments subject to a regularization based on the matrix information-entropy cutu13nips (). Once a kernel between two configurations has been defined, it is then possible to introduce a kernel distance
that we will use as the metric for representing and clustering structures from a database.
As discussed in Ref. de+16pccp (), the choice of a SOAP kernel as the definition of an environment similarity provides at the same time great generality – it can be seamlessly applied to both molecules and solids – and elbowroom for fine-tuning – ranging from setting an appropriate cutoff distance to circumscribe an environment to the introduction of an alchemical similarity kernel that translates the notion that different chemical species can behave similarly with respect to the properties of interest.
i.2 Mapping the Structural Landscape of a Database
The dissimilarity between the atomic configurations in a database contains a large amount of information on the structural relations between the database items. However, this information is not readily interpretable, as it is encoded as a matrix of numbers. Several methods are available to process dissimilarity information into a form that can be understood more intuitively. A first approach involves building a low-dimensional “map”, where each point corresponds to one of the structures in the database and where the (Euclidean) distances between points represents the information on the pairwise dissimilarity matrix.
Several methods have been proposed over the years to solve this dimensionality reduction problem, starting from principal component analysis pca-WOLD198737 () and the equivalent linear multi-dimensional scaling Kruskal1964 (), and proceeding to non-linear generalizations of the idea, such as ISOMAP tene+00science (), diffusion maps coif+05pnas (), kernel PCA scho+98nc (). In this manuscript, we will use sketchmap ceri+11pnas (); trib+12pnas (); ceri+13jctc (), a method in which one iteratively optimizes the objective function
that measures the mismatch of the dissimilarity between atomic configurations with the dissimilarity (typically just the Euclidean distance) between the corresponding low-dimensional projections . The procedure is very similar to multi-dimensional scaling, except for the appearance of the transformations and , which are non-linear sigmoid functions of the form:
The non-linear transformation focuses the optimization of Eq. 4 on the most significant distances (typically those of the order of ), and disregards local distortions (e.g. induced by thermal fluctuations or by incomplete convergence of a geometry optimization) and the relation between completely unrelated portions of configuration landscape. The maps that we report in this work will be labeled synthetically using the notation -A_B-a_b, where and denote the exponents used for the high-dimensional function , and denote the exponents for the low-dimensional function , and the threshold for the switching function. The choice of these parameters of the sigmoid functions are discussed in detail elsewhere ceri+13jctc (). In practice , , and have relatively small effect on the projection, and can be optimized and kept fixed for systems belonging to the same family. Since the structures we consider here are minimum-energy configurations, and there are no thermal fluctuations that should be filtered out, we set (so that at short range the algorithm will still try to represent distances faithfully) and set the long-range exponents to . The parameter is the one to which sketchmap is most sensitive, and needs to be tuned for each system separately. To automate the process of building sketchmaps of large amount of subsets of the database, we have used a simple heuristic procedure for determining the value of automatically. Following the prescriptions in Ref. ceri+13jctc (), we first compute the histogram of distances in the dissimilarity matrix of each molecular set, and detect the dissimilarity value () corresponding to the peak value of the histogram. We then set the value of to .
i.3 Hierarchical Clustering Representation
As we will demonstrate below, sketchmap provides a remarkably informative two-dimensional representation of structures in a data set, making it possible to identify groups of similar configurations, outliers, as well as to investigate structure-property relations. An alternative approach to navigate a set of structures based on the dissimilarity matrix is to use clustering algorithms, that identify groups of objects having similar properties to hint at the presence of recurring motifs underlying the behavior of the system.
A considerable number of clustering algorithms have been developed over the last few decades Jain:1999:DCR:331499.331504 (); clustering-rui (); aggarwal2013data (), including connectivity models WIDM:WIDM53 () (i.e. hierarchical clustering), centroid models kmeans1 (); kmeans2 (); kmeans3 () (i.e. k-means algorithm) and density based models dbscan (); optics (); rodr-laio14science ().
Clustering models based on connectivity information such as hierarchical (or agglomerative) clustering WIDM:WIDM53 () are particularly suited for this purpose and we will focus only on this type of clustering in this paper. Starting from each configurations as its own cluster, the hierarchical clustering algorithms iteratively aggregate clusters together based on some assessment of their similarity. Similarity between two individual structures can be obviously measured by their distance . The distance between two clusters, however, can be defined in many different ways. In our study, we will use in particular the RMS dissimilarity between the pair of members of the two clusters. The linkage distance between two clusters and is then defined as:
where and are the total number of configurations within each cluster. is the dissimilarity between the two configurations, as computed e.g. from the REMatch-SOAP kernel. The complexity of this type of clustering () is relatively cheaper than dimensionality reduction algorithms like sketchmap (). Both procedures can be greatly accelerated by selecting a subset of the configurations (e.g. by farthest point sampling, with the possibility of weighting based on density information ceri+13jctc ()) which is then subject to either dimensionality reduction or hierarchical clustering.
The results of a hierarchical clustering procedure can be represented in a “dendrogram” plot, that conveys visually the sequence of agglomerative clustering operations and the linkage distance at each step. The x-axis represents the index of the structures, sorted in a way so that at each stage the clusters being joint are adjacent. Each merge operation is represented by a line joining the two underlying clusters, with the position of the line representing the linkage distance for that pair, as defined by Eq. 6. In this kind of representation, at the bottom of the dendrogram, each structure can be thought of as an individual cluster containing only one item. Clusters are then merged iteratively, selecting at each step the pair of clusters that are closest to each other. This operation is repeated until all the clusters collapse into one single group that encompasses all the structures in the database, thus completing the dendrogram. To avoid overcrowding the bottom of the plot, one can hide the part that corresponds to very small linkage distances, while still graphically visualising the size of the clusters by drawing bars that encompass the associated structures. Since the “leaves” of this dendrogram correspond to individual configurations, it is possible to complement the dendrogram with color-coded bar plots that represent the value of different properties of each structure, thereby giving a clear picture of the relation between structural clustering and the different properties.
In order to understand the basic motifs of a particular cluster , it is very useful to select one of its structures that is as representative as possible of the entire subset. In case where stability estimates are available, such structure may be the lowest-energy structure in the cluster. For a definition that is based purely on conformational or configurational information, the most representative structure could be defined, as the item having the minimum mean square dissimilarity with respect to all other members of , i.e.
Representative structures can be defined at each level of the hierarchy, and can therefore be very useful in navigating the database, and understanding what are its most crucial structural features. The spread of the cluster around ,
can be used to quantify the range of structural landscape that is covered by the cluster.
Another important aspect of database analysis is ‘outlier detection’ outlier1 (); outlier2 (); outlier3 (); outlier4 (); outlier5 (); outlier6 (). An “outlier” configuration is defined as a configuration which is different from most of the configurations in the database. Outlier configurations are very important as they are likely to have unique structural motif in the whole database and are thus interesting for structure prediction applications. They also could represent chemical changes or indicate inconsistent configurations which are likely to be “errors” in the database.
In the following sections, we will present examples of how these different analyses can be applied to different subsets of structures taken from a database of amino acid and dipeptide conformers.
Ii Analysis of a database
This work is based on a first-principles derived structure/energy data set with conformers of twenty proteinogenic amino acids and dipeptides, as well as their interactions with a series of divalent cationsaminodb () (Ca, Ba, Sr, Cd, Pb, Hg). The potential-energy surfaces (PES) of 280 systems were explored in a wide relative energy range of up to 4 eV (390 kJ/mol), summing up to an overall of 45,892 stationary points on the respective potential-energy surfaces db-paper (). The underlying energetics were calculated by applying density-functional theory (DFT) in the generalized gradient approximation corrected for long-range van der Waals interactionscpc180_2175 (); prl77_3865 (); prl102_73005 () (PBE+vdW). A number of theory-theory and theory-experiment comparisons have shown the applicability of the method to amino acid and peptide systems prl106_118102 (); cej19_11224 (); db-paper (); C4CP05541A (); C4CP05216A (); doi:10.1021/jp412055r (); 0953-8984-27-49-493002 (). Such dominantly manually curated data is a great starting point for studying materials and molecules across chemical space ropo_scirep_2016 (), yet automated techniques are needed to extract unbiased and hypothesis-free trends and to check consistency of the underlying data.
In this study we focus on the amino acid lysine (in short Lys) and investigate basic structural motifs of three forms, see Fig. 1. Furthermore, the machine learning techniques introduced in this work are used to detect the impact of perturbations (here Ca cations) on the structural properties of the unperturbed systems. Finally, we demonstrate how the approach can also be applied to discover inconsistencies and outliers in the database. Hierarchical classifications and sketchmap projections for all the proteinogenic amino acids in the database are given in the Supporting Information.
ii.1 Finding the Dominant Features of a Structual Landscape
ii.1.1 Lysine Dipeptide
We take as our first example a subset of the database containing 2080 conformers of lysine dipeptide. As discussed in the previous section, we start by constructing the (dis)similarity matrix using the SOAP-REMatch kernel. In Fig. 2 the dendrogram plot as well as sketchmaps have been shown along with five properties, energy and four dihedral angles, using the same color scales in both the sketchmap and dendrogram representations. In the sketchmap each circular ‘disk’ represents a conformer. Whereas in the case of the dendrogram plot, structures are represented by vertical lines at the bottom of the plot. The strong correlation between energy and conformational parameters on one side, and clustering and position on the map on the other, testifies how the the REMatch-SOAP kernel induces a meaningful classification of the structures in this dataset.
While both clustering and sketchmap show clearly that the dataset is composed of groups of structurally-related conformers, the agnostic nature of the underlying metric does not disclose immediately the structural features that most transparently differentiate between different clusters. Comparing the representative structures from the main clusters allowed us to quickly identify candidate structural motifs that could be used to rationalize the layout of the conformational landscape. By color-coding the dendrogram and the sketchmaps according to these indicators one can readily highlight the key correlations. When considering existing literature on the stability of oligopeptides, the two structural parameters that are most often considered as the key coordinates to navigate the conformational landscape are the Ramachandran dihedral angles \straightphi and \textpsi, that determine the structure of the backbone around the side chain bearing C\textalpha atom ramachandran-plot () under the assumption of peptide bonds being solely in trans conformation. While fine-grained clusters are clearly homogeneous with respect to the \straightphi and \textpsi angles, it is clear that for the present systems they are not those that determine the clear-cut branching at the top of the dendrogram. An analysis of the representative structures for the two main clusters (a) and (b) shows that the two molecules differ by the isomerization of the N-terminal peptide bond. Further splitting of these two clusters, i.e. (a) into clusters (d) and (e), and (b) into (f) and (g), depends on the isomerization of the C-terminal peptide bond. We can confirm this attribution of the main features of the dataset by color-coding the map and the dendrogram following the dihedral angles \textomega and \textomega. The four main clusters are largely homogeneous with respect to peptide bond isomerization, and are then further subdivided based on \straightphi and \textpsi. This observation deserves some further comment. Peptide bonds in naturally-occurring proteins are believed to almost exclusively exist in trans conformation with the exception of prolyl peptide bonds where a smaller energy difference to trans increases the chance for cis onformersfischer_chemical_2000 (); dugave_cistrans_2003 (). This view is supported by the analysis of protein structures deposited in the protein databases where cis conformations are found for about 5% of the prolyl peptide bonds, but less than 0.1% for the othersweiss_peptide_1998 (). X-ray crystallographic structure represent however merely frozen snapshots of structural dynamics. The ab initio structure search protocol, instead, does consider the peptide bond torsions as variable and intentionally allowed simulations to overcome the isomerization barrier. Consequently, the dataset contains representatives of all four combinations of cis and trans conformers. Since these transitions are strongly bimodal, and reflect in significant changes of the favorable side chain conformations, they constitute the most significant feature to classify the conformers. As expected, the most stable conformers are largely in a trans-trans conformation. However, the large parts of conformational space of that is occupied by conformers with 1 or 2 cis peptide bonds suggests that cis isomers might play a role in the dynamics of peptides and proteins. If the analysis had been performed solely under the assumption of the importance of the Ramachandran \straightphi and \textpsi dihedrals, it would have missed one of the main features of the structural landscape that is critical to characterize the relation between structure and energetics. One could then proceed further with the analysis, focusing for instance on small clusters containing low-energy structures such as that represented by the conformer (c). All the structure in this group are trans-trans isomers, that in addition have degrees and degrees, allowing for the formation of a H-bond between the side chain and , and a favorable arrangement of the donating a H-bond to the carbonyl as shown in Fig. 2. Having access to the combined information on energetics, and on the grouping of structures with similar geometry makes it easier to rationalize the energy ordering of the structures, without having to separately juxtapose all of the low-lying conformers but focusing on a few representative structures.
ii.1.2 Protonated Lysine Dipeptide
As the second example we considered a dataset containing 897 conformers of gas phase protonated lysine dipeptide. We follow the same steps as described in the previous example in order to find the most basic structural motifs for this system. Fig. 3 shows the dendrogram, the sketchmap and a few color coded properties for this system to show their correlation with the classification. The most prominent feature for this molecule, which is evident in both the dendrogram and the sketch maps, is the presence of a group of outliers, that are clearly separated from the bulk of the conformers. Inspection of the cluster centroid (g) clarifies the structural basis of this separation. Whereas in most of the structures the excess charge lies on the lysine side chain as a NH group, conformers in this cluster experienced a proton transfer event, with the excess proton attached to one of the carbonyl oxygen , stabilized by H-bonding to . This is a result of the database generation where ab initio replica-exchange molecular dynamics including high T trajectories where used for structure sampling during which protons can eventually transfer.
Moving on to the main cluster of structures, we can see that similar to our previous example of the neutral dipeptide and again due to the unbiased sampling protocol and the high energy range the peptide bond angles are again more important than Ramachandran’s dihedrals. Conformers (a) and (b) are the representative structure for groups having cis and trans \textomega peptide bonds respectively. Group (a) is further split based on the cis/trans state of \textomega into the clusters represented by structures (d) and (e).
The presence of a charged side chain leads to stronger H-bonds and to less clear-cut grouping based on peptide-bond isomerism than for the neutral dipeptide. This is apparent for instance from the presence of subclusters such as that represented by conformer (f), in which the bent side chain leads to the formation of two H-bonds between NH group and the carbonyl oxygens. H-bonds also dominate the partitioning of cluster (b), that is split into two groups – one of which is still best represented by the same conformer, and one that is epitomised by (c). Once again, inspection of these structural representatives reveals the organising principle behind the classification: (c)-like structures have an extended side chain, and are dominated by interactions among the peptide bond moieties, whereas (b)-like structures have a well-formed H-bond between the side chain and one of the two backbone O atoms. This structural pattern can be emphasized by color-coding conformers based on the parameter : A small O-N distance indicates bending of the side chain and the formation of a H-bond between O and N. As it is clear from the sketchmap representation, there is a very strong correlation between the bending of the charged side chain and the energy of a conformer, with all of the structures within 0.5 eV of the ground state feature this sidechain to backbone H-bonds.
A scenario with such a clear-cut partitioning is less likely to happen with oligopeptides in a polar solvent like water, where intramolecular H-bonds that introduce strain compete with H-bonds with the surrounding water molecules, without the need of bending the structure. The analysis techniques we introduce in this work would be ideally suited to rationalize the changes in the (free) energetics of biological molecules when moving from the gas phase to (micro)solvated environments or to organic/inorganic interfaces.
ii.1.3 Uncapped Lysine
Our third example is a dataset containing 733 conformers of the un-capped lysine molecule in the gas phase. We follow the same steps as described in the previous examples to construct the dendrogram shown in Fig. 4. The map has a simple structure, with few well-separated groups. Being a smaller molecule with fewer degrees of freedom, the Ramachandran angles are not defined. Still, the dihedral angles in the vicinity of the atom display local structural correlation but once again they are not the main organizing factor that can rationalize the clustering. By juxtaposing representative conformers from the main clusters we could identify a better order parameter, that correlates strongly with H-bond patterns within the molecule. Namely, the distance () between the H atom in the group of the carboxy function and the atom in the backbone () discriminates well between structures based on H-bonding patternsropo_scirep_2016 () of type I between (e.g. conformer (b)) and of type II with a H-bond (e.g. conformer (a)). It can be seen from both the dendrogram and the sketchmaps that one could identify several further subgroups based on particular values of , representing specific orientations. Conformers (c) and (d) represent small groups of conformers having specific relative orientation between the OH and groups. Conformer (e) is representative of a small outlier group with a well-defined bend of the side chain, leading to the formation of a further H-bond between the atom in the amino acid moiety and , in the side chain. The lysine side chain is very flexible, and the distance between and only plays a role in defining the fine-grained structure of the dataset, but is minimally correlated with the most prominent features.
While it appears that even in this case we could identify the basic structural motifs that characterize the conformational landscape of this system, the correlation with energy is very poor. There are several instances in both the dendrogram and the sketchmap where two conformers that are detected as structurally very similar while they display very different stability. Understanding whether this inconsistency signals a problem with our analysis brings us to the topic of outlier detection and consistency checks, that we will discuss in details in Section II.3.
ii.2 Understanding the Impact of Perturbations on conformational Space
Having elucidated the essential structural motifs that underlie the organization of a set of molecular conformers, one could also wonder how changes in the thermodynamic conditions, or another external perturbation such as solvation, the addition or subtraction of an electron Dejcp11 () or that of an atom ideh_mg (); ghazi_Na (); pascalB80 (), modify the conformations of the molecule and their stability. The databasedb-paper (); aminodb () that we are using here as an example contains, in addition to bare oligopeptides, sets of locally-stable conformers in the presence of cations of six different species, namely \ceCa^2+, \ceBa^2+, \ceSr^2+, \ceCd^2+, \cePb^2+ and \ceHg^2+. We take the example of \ceCa^2+ to describe how one can characterize its impact on the conformational space of the three molecular systems that we have discussed in our previous examples. We start by calculating the dissimilarity of all the conformers containing cations with their pure counterpart. In order to make the comparison on the same footings, we did not include the location of the cation in defining the SOAP kernels, so that the presence of \ceCa^2+ only enters by distorting molecular geometries and/or altering their relative stability. Using this information, we then projected the cation-containing dataset on the top of the sketchmap of structures for the bare molecule. This is done using sketchmap out-of-sample embedding, and we refer our reader to see the relevant literature ceri+11pnas (); trib+12pnas (); ceri+13jctc () for more details about the method. In Fig. 5 we show the resulting projection, colored according to the stability of the conformers, on top of the sketchmap of the pure molecule shown in grey color as a reference. A close proximity of projected conformers with a pure conformer signifies their structural similarity. Segregation of the projected conformers with the cation in some area of the reference sketchmap, represents the structural bias introduced by the strong electrostatic interaction with \ceCa^2+.
In the case of neutral lysine dipeptide (Fig. 5-a), the presence of the Ca ion induces overall relatively small distortions of the stable conformers, that get pushed towards the outer region of the map but are still clearly related to the locally stable structures for the bare molecule. Energies are dramatically changed, with the most stable cluster in the original map begin completely absent in the presence of the cation. These observations highlight the importance of sampling high-energy conformers during high-throughput structure searches, since the relative stability can be modulated strongly by external perturbations. In particular, cis conformers become energetically more competitive and are topologically closer to the global minima. In the case of protonated lysine dipeptide (Fig. 5-b), the same analysis shows an even clearer pattern. All the conformers with \ceCa^2+ ions are projected in the lower part of the sketchmap, that correspond to conformers with an extended side chain (see Fig. 3). The \ceCa^2+ ion preferably binds to the peptide atoms, and the electrostatic repulsion with the protonated lysine residue strongly favors extended conformers, contrary to what observed in the case of the bare molecule. Finally, one sees that for molecular lysine (Fig. 5-c) the addition of Ca leads to conformers with very different structural motifs from those seen with the bare molecule, which is apparent in the sketchmap projection being concentrated far away from the unperturbed conformers. In fact, inspection of the structures shows that in most cases Ca triggers the transition to the zwitterionic form, with the cation coupled to the carboxylate group, and the side chain protonated NH extending as far as possible away from it. In analogy with what observed for Lennard-Jones clusters ceri+13jctc () and solvated polypeptide segments arde+15jctc (), sketchmaps proved to be a powerful tool to analyze the response of the system to external perturbations and changes in the boundary conditions, and – in this specific example – do draw connections between different subsets of a high-throughput molecular database.
ii.3 Identifying Outliers and Checking for Consistency
The tools we introduced in this work are useful to address other important issues in data-driven science, such as outlier detection and consistency checks. We have already discussed the importance of detecting groups of “outlier” structures that are very different from the bulk of the dataset. These unusual items often signal the occurrence of unexpected effects that go beyond the original goal of the database construction effort. In the case of protonated lysine dipeptide, looking for outliers allowed us to reveal the presence of conformers with different chemical connectivity, or of strong H-bonds between the backbone and the charged side chain. Similar observations can also be made in the case of the bare lysine molecule (Fig. 4). Also in this case, one can observe a branch at the topmost level of the dendrogram, containing only two conformers. These are the only two cases where a H-bond is formed between the of the side chain and the atom of the group in the backbone. In the sketchmap, these two conformers are projected on the top, clearly isolated from rest of the groups, and bear the most resemblance to the zwitterionic conformers that are stabilized in the presence of a divalent cation. Obviously, the definition of a group of “outliers” can be more nuanced, and refer to small groups of structures appearing at deeper levels in the hierarchy. Overall, the possibility of clustering together the structures in a large dataset and inspect a few representative conformers, rather than hundreds or thousands, greatly facilitates the task of identifying trends and spotting interesting or unexpected structures.
Outliers can signal interesting or important trends, but can also be a red flag for the presence of errors. The importance of database integrity has long been recognized by computer scientists dataint1 (); dataint2 (); dataint3 (); dataint4 (), and several tools are available to monitor and correct inconsistencies from the technical point of view, in terms of reliability of storing and retrieving dataoutlier1 (); outlier2 (); outlier3 (); outlier4 (); outlier5 (); outlier6 (). The issue is also crucial when it comes to the maintenance of automatically-generated databases, and to repositories that store data of heterogeneous provenanceaiida (); clean_en (); esp (); oqmd (); pauling (); matgenome1 (); matgenome2 (). In these cases, problems have generally little to do with the integrity of the storage, but rather with the consistency of the simulation details of different sets of calculations. Inconsistencies should here manifest themselves in the presence of structures that are geometrically very similar, but are associated to very different values of particular properties.
For example the lysine molecule dataset shows signs of this kind of issues, with energies that vary wildly within clusters that are very homogeneous in structure. This problem can be seen from the maps, i.e. when comparing the energy-colored sketchmap in Fig. 4 to the respective maps for the other systems. However, a more robust and easy-to-automate approach to identify structure/property inconsistencies starts from the hierarchical clusters, and compares the structural variability within each cluster (Eq. 8) with the variance of a given property, in this case energy, . Looking, for example, at a glassy energy landscapeDe2011 (), one can observe configurations that are very different from a structural point of view, but have similar energy, giving rise to clusters with large and small . The data points in Fig. 6 each represent individual clusters of lysine dipeptide and uncapped lysine, respectively, and illustrate their variation in energy and structure. In the case of lysine dipeptide (Fig. 6a) one sees a clear correlation between the structural and energetical variation of the clusters. The two quantities and are not necessarily strongly correlated, but in general clusters that contain very similar structures also have a low spread in energy. For uncapped lysine (Fig. 6b), however, one can identify a group of points (which we manually highlighted in orange for clarity) that has a distinctively different behavior, with converging to a constant value other than zero as decreases. This kind of feature indicates that the metric that is used to classify structures cannot detect one specific effect that has a dramatic impact on energetics, signaling either a failure of the metric or, as in this case, an inconsistency in the generated data. Further investigation of the lysine molecule dataset revealed that a subset of structures that had been generated at a lower level of theory in the initial stages of the structure-search procedure made their way by mistake into the final dataset. Using this measure of cluster homogeneity on all systems of the amino acid database (see Supporting Information) revealed similar problems also for other molecules, for example Cys, Glu, and Arg. Thanks to this analysis we will be able to identify and rectify mistakes in all the affected datasets and subsequently update the on-line repositoryaminodb ().
The increasing use of high-throughput computational screening of materials and molecules, and the compilations of more-or-less curated databases of the resulting structures and properties, is making more and more urgent to adapt “big data” techniques to the problems that are specific to this field. In this work we have demonstrated how a toolbox of algorithms ranging from hierarchical clustering to non-linear dimensionality reduction can be used to navigate molecular databases, taking as a paradigmatic example some subsets of a database of oligopeptide structures in the gas phase. The software that was used to compute similarity data between molecules, as well as to generate dendrograms and sketch-maps, are open-source and available for downloadcosmogit (); libatoms ().
We find that the use of REMatch-SOAP, a general and unbiased metric to compare different structures based on the combination of pair-wise similarity between molecular environments, makes these techniques particularly insightful. Rather than simply reflecting preconceived notions of what would be the key structural parameters to differentiate molecular conformers, this metric reveals for instance the importance of peptide bond isomerization in describing the high-energy portion of conformational space of oligopeptides, the possibility of changes in chemical connectivity in the course of the ab initio structural search, and the interplay between hydrogen-bonding, backbone dihedrals an electrostatic interactions. Sketchmaps and hierarchical clustering proved to be complementary tools, with representative structures from the main clusters providing an easy way to compare visually different groups of conformers, and the low-dimensional map providing a quick, intuitive tool to verify hypotheses and visualize structure-property correlations.
Assumption-free first-principles molecular-structure search for data generation in combination with dimensionality reduction and clustering for data analysis provide a powerful tool box for studying structure formation trends. We could highlight the presence of large portions of configurational space that consist of cis isomers of the peptide bond. Albeit energetically unfavorable, these conformers may play an important role in the dynamics of polypeptides. By comparing isolated molecules with their complexes with \ceCa^2+, we can also reveal how a strong electrostatic perturbation modifies the energetic landscape of a small molecule – be it by shifting the stability of different conformers, or triggering the formation of new structures that are absent in the absence of a cation. Furthermore, we also demonstrate the importance of automated analysis techniques in assessing the integrity and internal consistence of a database, by successfully identifying a subset of structures associated with ill-converged energetics.
All of the techniques we discussed should be readily extendable to heterogeneous databases of molecules and solids, where we expect that the possibility of defining an alchemical kernel within the REMatch-SOAP metric will make it possible to tune the relative weight of composition and structure in determining the notion of similarity. By simplifying the analysis and the interpretation of computational datasets containing thousands or millions of hypothetical compounds, these methods will be crucial to unleash the full potential of computational materials design.
iv.1 Availability of data and materials
The oligopeptide database used as an example in this paper is already available online aminodb (). All the data required to generate the figures in this paper and same analysis data for other oligopeptides in the database are provided in the supplementary information. Similarity matrix data is not included in SI due to size limitations, but is freely available from from our group repository cosmogit ().
iv.2 Competing interests
We confirm that none of the authors have any competing interests in the manuscript.
S.D. and M.C. would like to acknowledge support from the NCCR MARVEL. M.C., T.I. and C.B. would like to acknowledge funding from the MPG-EPFL center for molecular nanoscience.
iv.4 Authors’ contributions
M.C. and S.D. designed the calculations and developed the methods. F.M. and S.D. performed calculations and analysis, and prepared the materials for the manuscript. C.B. and T.I. provided insights on the implications of the analysis of the oligopeptide database. All the authors contributed to the writing of the manuscript.
S.D. would like to thank Czuee Morey (University of Geneva) for insightful discussion. C.B. thanks Matti Ropo (Tampere University of Technology) Volker Blum (Duke University) and Matthias Scheffler (Fritz Haber Institute) for support and discussion.
- (1) Pizzi G, Cepellotti A, Sabatini R, Marzari N, Kozinsky B. AiiDA: automated interactive infrastructure and database for computational science. Computational Materials Science. 2016;111(1):218–230.
- (2) Hachmann J, Olivares-Amaya R, Atahan-Evrenk S, Amador-Bedolla C, Sánchez-Carrera RS, Gold-Parker A, et al. The harvard clean energy project: Large-scale computational screening and design of organic photovoltaics on the world community grid. Journal of Physical Chemistry Letters. 2011;2(17):2241–2251.
- (3) Ortiz C, Eriksson O, Klintenberg M. Data mining and accelerated electronic structure theory as a tool in the search for new functional materials. Computational Materials Science. 2009;44(4):1042–1049.
- (4) Saal JE, Kirklin S, Aykol M, Meredig B, Wolverton C. Materials Design and Discovery with High-Throughput Density Functional Theory: The Open Quantum Materials Database (OQMD). JOM. 2013;65(11):1501–1509.
- (5) Villars P, Berndt M, Brandenburg K, Cenzual K, Daams J, Hulliger F, et al. The Pauling File, Binaries Edition. Journal of Alloys and Compounds. 2004;367(1-2):293–297.
- (6) Jain A, Ong SP, Hautier G, Chen W, Richards WD, Dacek S, et al. Commentary: The materials project: A materials genome approach to accelerating materials innovation. APL Materials. 2013;1(1):011002.
- (7) White A. The Materials Genome Initiative: One year on. MRS Bulletin. 2012;37(08):715–716.
- (8) Rupp M, Tkatchenko A, Müller KR, von Lilienfeld OA. Fast and Accurate Modeling of Molecular Atomization Energies with Machine Learning. Physical Review Letters. 2012;108(5):058301.
- (9) Ghiringhelli LM, Vybiral J, Levchenko SV, Draxl C, Scheffler M. Big data of materials science: Critical role of the descriptor. Physical Review Letters. 2015;114(10):105503.
- (10) Huan TD, Mannodi-Kanakkithodi A, Ramprasad R. Accelerated materials property predictions and design using motif-based fingerprints. Physical Review B - Condensed Matter and Materials Physics. 2015;92(1):14106.
- (11) Botu V, Ramprasad R. Learning scheme to predict atomic forces and accelerate materials simulations. Physical Review B - Condensed Matter and Materials Physics. 2015;92(9):94306.
- (12) Kusne A, Gao T, Mehta A, Ke L, Cuong Nguyen M, Ho KM, et al. On-the-fly machine-learning for high-throughput experiments: search for rare-earth-free permanent magnets. Scientific Reports. 2014;4:6367.
- (13) Ramakrishnan R, Dral PO, Rupp M, von Lilienfeld OA. Quantum chemistry structures and properties of 134 kilo molecules. Scientific data. 2014;1:140022.
- (14) Arsenault LF, Lopez-Bezanilla A, Von Lilienfeld OA, Millis AJ. Machine learning for many-body physics: The case of the Anderson impurity model. Physical Review B - Condensed Matter and Materials Physics. 2014;90(15):155136.
- (15) Ropo M, Schneider M, Baldauf C, Blum V. First-principles data set of 45,892 isolated and cation-coordinated conformers of 20 proteinogenic amino acids. Scientific Data. 2016;3:160009.
- (16) Rodriguez A, Laio A. Clustering by fast search and find of density peaks. Science. 2014;344(6191):1492–1496.
- (17) Xu R, Wunsch D. Survey of clustering algorithms. IEEE Transactions on Neural Networks. 2005;16(3):645–678.
- (18) Yu G, Chen J, Zhu L. Data mining techniques for materials informatics: Datasets preparing and applications. In: 2009 2nd International Symposium on Knowledge Acquisition and Modeling, KAM 2009. vol. 2; 2009. p. 189–192.
- (19) Isayev O, Fourches D, Muratov EN, Oses C, Rasch K, Tropsha A, et al. Materials cartography: Representing and mining materials space using structural and electronic fingerprints. Chemistry of Materials. 2015;27(3):735–743.
- (20) Balachandran PV, Theiler J, Rondinelli JM, Lookman T. Materials Prediction via Classification Learning. Scientific Reports. 2015;5:13285.
- (21) Ferguson AL, Panagiotopoulos AZ, Debenedetti PG, Kevrekidis IG. Systematic determination of order parameters for chain dynamics using diffusion maps. Proceedings of the National Academy of Sciences of the United States of America. 2010;107(31):13597–602.
- (22) Ceriotti M, Tribello GA, Parrinello M. From the Cover: Simplifying the representation of complex free-energy landscapes using sketch-map. Proceedings of the National Academy of Sciences. 2011;108(32):13023–13028.
- (23) Tribello Ga, Ceriotti M, Parrinello M. Using Sketch-Map Coordinates to Analyze and Bias Molecular Dynamics Simulations. Proceedings of the National Academy of Sciences. 2012;109(14):5196–5201.
- (24) Ceriotti M, Tribello GA, Parrinello M. Demonstrating the transferability and the descriptive power of sketch-map. Journal of Chemical Theory and Computation. 2013;9(3):1521–1532.
- (25) Rohrdanz MA, Zheng W, Clementi C. Discovering Mountain Passes via Torchlight: Methods for the Definition of Reaction Coordinates and Pathways in Complex Macromolecular Reactions. Annual Review of Physical Chemistry. 2013;64(1):295–316.
- (26) De S, Bartók AP, Csányi G, Ceriotti M. Comparing molecules and solids across structural and alchemical space. Phys Chem Chem Phys. 2016;18(20):13754.
- (27) Ropo M, Baldauf C, Blum V. Berlin ab initio amino acid DB; 2016.
- (28) Pietrucci F, Andreoni W. Graph theory meets ab initio molecular dynamics: Atomic structures and transformations at the nanoscale. Physical Review Letters. 2011;107(8):85504.
- (29) Szlachta WJ, Bartók AP, Csányi G. Accuracy and transferability of Gaussian approximation potential models for tungsten. Physical Review B - Condensed Matter and Materials Physics. 2014;90(10):104108.
- (30) Lopez-Bezanilla A, Von Lilienfeld OA. Modeling electronic quantum transport with machine learning. Physical Review B - Condensed Matter and Materials Physics. 2014;89(23):235411.
- (31) Pilania G, Wang C, Jiang X, Rajasekaran S, Ramprasad R. Accelerating materials property predictions using machine learning. Scientific Reports. 2013;3:2810.
- (32) Bartók AP, Gillan MJ, Manby FR, Csányi G. Machine-learning approach for one- and two-body corrections to density functional theory: Applications to molecular and condensed water. Physical Review B - Condensed Matter and Materials Physics. 2013;88(5):054104.
- (33) Rupp M, Proschak E, Schneider G. Kernel approach to molecular similarity based on iterative graph similarity. Journal of Chemical Information and Modeling. 2007;47(6):2280–2286.
- (34) Hirn M, Poilvert N, Mallat S. Quantum Energy Regression using Scattering Transforms. arXiv preprint arXiv:150202077. 2015;.
- (35) Montavon G, Rupp M, Gobre V, Vazquez-Mayagoitia A, Hansen K, Tkatchenko A, et al. Machine learning of molecular electronic properties in chemical compound space. New Journal of Physics. 2013;15(9):95003.
- (36) Snyder JC, Rupp M, Hansen K, Müller KR, Burke K. Finding density functionals with machine learning. Physical Review Letters. 2012;108(25):253002.
- (37) Ghasemi SA, Hofstetter A, Saha S, Goedecker S. Interatomic potentials for ionic systems with density functional accuracy based on charge densities obtained by a neural network. Physical Review B. 2015;92(4):045131.
- (38) Von Lilienfeld OA. First principles view on chemical compound space: Gaining rigorous atomistic control of molecular properties. International Journal of Quantum Chemistry. 2013;113(12):1676–1689.
- (39) Hansen K, Biegler F, Ramakrishnan R, Pronobis W, Von Lilienfeld OA, Müller KR, et al. Machine learning predictions of molecular properties: Accurate many-body potentials and nonlocality in chemical space. Journal of Physical Chemistry Letters. 2015;6(12):2326–2331.
- (40) Zhu L, Amsler M, Fuhrer T, Schaefer B, Faraji S, Rostami S, et al. A fingerprint based metric for measuring similarities of crystalline structures. The Journal of Chemical Physics. 2016;144(3):034203.
- (41) Cuturi M. Sinkhorn Distances: Lightspeed Computation of Optimal Transport. In: Burges CJC, Bottou L, Welling M, Ghahramani Z, Weinberger KQ, editors. Advances in Neural Information Processing Systems 26. Curran Associates, Inc.; 2013. p. 2292–2300.
- (42) Wold S, Esbensen K, Geladi P. Principal component analysis. Chemometrics and Intelligent Laboratory Systems. 1987;2(1):37–52.
- (43) Kruskal JB. Nonmetric multidimensional scaling: A numerical method. Psychometrika. 1964;29(2):115–129.
- (44) Tenenbaum JB, de Silva V, Langford JC. A global geometric framework for nonlinear dimensionality reduction. Science (New York, NY). 2000;290(5500):2319–23.
- (45) Coifman RR, Lafon S, Lee aB, Maggioni M, Nadler B, Warner F, et al. Geometric diffusions as a tool for harmonic analysis and structure definition of data: diffusion maps. Proceedings of the National Academy of Sciences of the United States of America. 2005;102(21):7426–31.
- (46) Schölkopf B, Smola A, Müller KR. Nonlinear Component Analysis as a Kernel Eigenvalue Problem. Neural Computation. 1998;10(5):1299–1319.
- (47) Jain, A K, Murty, M P, Flynn, P J. Data clustering: a review. ACM Computing Surveys. 1999;31(3):264–323.
- (48) Aggarwal CC, Reddy CK. Data Clustering: Algorithms and Applications. CRC Press; 2013.
- (49) Murtagh F, Contreras P. Algorithms for hierarchical clustering: an overview. Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery. 2012;2(1):86–97.
- (50) Huang Z. Extensions to the k-means algorithm for clustering large data sets with categorical values. Data Mining and Knowledge Discovery. 1998;2(3):283–304.
- (51) Jing L, Ng MK, Huang JZ. An entropy weighting k-means algorithm for subspace clustering of high-dimensional sparse data. IEEE Transactions on Knowledge and Data Engineering. 2007;19(8):1026–1041.
- (52) Su MC, Chou CH. A modified version of the K-means algorithm with a distance based on cluster symmetry. IEEE Transactions on Pattern Analysis and Machine Intelligence. 2001;23(6):674–680.
- (53) Ester M, Kriegel HP, Sander J, Xu X. A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise. In: Proceedings of the 2nd International Conference on Knowledge Discovery and Data Mining. AAAI Press; 1996. p. 226–231.
- (54) Ankerst M, Breunig MM, Kriegel HP, Sander J. Optics: Ordering points to identify the clustering structure. In: ACM Sigmod Record. ACM Press; 1999. p. 49–60.
- (55) Zhao X, Liang J, Cao F. A simple and effective outlier detection algorithm for categorical data. International Journal of Machine Learning and Cybernetics. 2014;5(3):469–477.
- (56) Yamanishi K, Takeuchi JI, Williams G, Milne P. On-line unsupervised outlier detection using finite mixtures with discounting learning algorithms. Data Mining and Knowledge Discovery. 2004;8(3):275–300.
- (57) Petrovskiy MI. Outlier detection algorithms in data mining systems. Programming and Computer Software. 2003;29(4):228–237.
- (58) Angiulli, Fabrizio and Pizzuti C. In: Elomaa T, Mannila H, Toivonen H, editors. Fast Outlier Detection in High Dimensional Spaces. vol. 2431. Berlin, Heidelberg: Springer Berlin Heidelberg; 2002. p. 15–27.
- (59) Breunig MM, Kriegel HP, Ng RT, Sander J. LOF: identifying density-based local outliers. ACM SIGMOD Record. 2000;29(2):93–104.
- (60) Aggarwal CC, Yu PS, Aggarwal CC, Yu PS. Outlier detection for high dimensional data. Proceedings of the 2001 ACM SIGMOD international conference on Management of data - SIGMOD ’01. 2001;30(2):37–46.
- (61) Blum V, Gehrke R, Hanke F, Havu P, Havu V, Ren X, et al. Ab initio molecular simulations with numeric atom-centered orbitals. Computer Physics Communications. 2009;180(11):2175–2196.
- (62) Perdew JPJ, Burke K, Ernzerhof M, of Physics D, Quantum Theory Group Tulane University NOLJ. Generalized Gradient Approximation Made Simple. Physical Review Letters. 1996;77(18):3865–3868.
- (63) Tkatchenko A, Scheffler M. Accurate Molecular Van Der Waals Interactions from Ground-State Electron Density and Free-Atom Reference Data. Physical Review Letters. 2009;102(7):073005.
- (64) Tkatchenko A, Rossi M, Blum V, Ireta J, Scheffler M. Unraveling the stability of polypeptide helices: Critical role of van der Waals interactions. Physical Review Letters. 2011;106(11):118102.
- (65) Baldauf C, Pagel K, Warnke S, Von Helden G, Koksch B, Blum V, et al. How cations change peptide structure. Chemistry - A European Journal. 2013;19(34):11224–11234.
- (66) Schubert F, Rossi M, Baldauf C, Pagel K, Warnke S, von Helden G, et al. Exploring the conformational preferences of 20-residue peptides in isolation: Ac-Ala19-Lys + H(+)vs. Ac-Lys-Ala19+ H(+) and the current reach of DFT. Phys Chem Chem Phys. 2015;17(11):7373–7385.
- (67) Schubert F, Pagel K, Rossi M, Warnke S, Salwiczek M, Koksch B, et al. Native like helices in a specially designed peptide in the gas phase. Phys Chem Chem Phys. 2015;17(7):5376–5385.
- (68) Rossi M, Chutia S, Scheffler M, Blum V. Validation Challenge of Density-Functional Theory for Peptides-Example of Ac-Phe-Ala5-LysH(+). The journal of physical chemistry A. 2014;118(35):7349–59.
- (69) Baldauf C, Rossi M. Going clean: structure and dynamics of peptides in the gas phase and paths to solvation. Journal of physics Condensed matter : an Institute of Physics journal. 2015;27(49):493002.
- (70) Ropo M, Blum V, Baldauf C. Trends for isolated amino acids and dipeptides: Conformation, divalent ion binding, and remarkable similarity of binding to calcium and lead. arXiv:160602151 [physics, q-bio]. 2016;.
- (71) Ramachandran GN, Ramakrishnan C, Sasisekharan V. Stereochemistry of polypeptide chain configurations. Journal of molecular biology. 1963;7(1):95–99.
- (72) Fischer G. Chemical aspects of peptide bond isomerisation. Chemical Society Reviews. 2000;29(2):119–127.
- (73) Dugave C, Demange L. Cis-trans isomerization of organic molecules and biomolecules: Implications and applications. Chemical Reviews. 2003;103(7):2475–2532.
- (74) Weiss MS, Jabs A, Hilgenfeld R. Peptide bonds revisited. Nature structural biology. 1998;5(8):676.
- (75) De S, Ghasemi SA, Willand A, Genovese L, Kanhere D, Goedecker S. The effect of ionization on the global minima of small and medium sized silicon and magnesium clusters. Journal of Chemical Physics. 2011;134(12):124302.
- (76) Heidari I, De S, Ghazi SM, Goedecker S, Kanhere DG. Growth and structural properties of MgN (N = 10-56) clusters: Density functional theory study. Journal of Physical Chemistry A. 2011;115(44):12307–12314.
- (77) Ghazi SM, De S, Kanhere DG, Goedecker S. Density functional investigations on structural and electronic properties of anionic and neutral sodium clusters Na N (N = 40-147): comparison with the experimental photoelectron spectra. Journal of Physics: Condensed Matter. 2011;23(40):405303.
- (78) Pochet P, Genovese L, De S, Goedecker S, Caliste D, Ghasemi SA, et al. Low-energy boron fullerenes: Role of disorder and potential synthesis pathways. Physical Review B - Condensed Matter and Materials Physics. 2011;83(8):81403.
- (79) Ardevol A, Tribello GA, Ceriotti M, Parrinello M. Probing the unfolded configurations of a -hairpin using sketch-map. Journal of Chemical Theory and Computation. 2015;11(3):1086–1093.
- (80) Baškarada S, Koronios A. A Critical Success Factor Framework for Information Quality Management. Information Systems Management. 2014;31(4):276–295.
- (81) Van Den Broeck J, Cunningham SA, Eeckels R, Herbst K. Data cleaning: Detecting, diagnosing, and editing data abnormalities. PLoS Medicine. 2005;2(10):0966–0970.
- (82) Gevorgyan A, Poolman MG, Fell DA. Detection of stoichiometric inconsistencies in biomolecular models. Bioinformatics. 2008;24(19):2245–2251.
- (83) Ferretti L, Colajanni M, Marchetti M. Distributed, concurrent, and independent access to encrypted cloud databases. IEEE Transactions on Parallel and Distributed Systems. 2014;25(2):437–446.
- (84) De S, Willand A, Amsler M, Pochet P, Genovese L, Oedecker S. Energy landscape of fullerene materials: A comparison of boron to boron nitride and carbon. Physical Review Letters. 2011;106(22):225502.
- (85) Code repositories from the Laboratory of Computational Science and Modelling at EPFL; 2014. http://epfl-cosmo.github.io/.
- (86) Libatoms; 2014. http://www.libatoms.org/.