Article pubs.acs.org/jcim
aRMSD: A Comprehensive Tool for Structural Analysis Arne Wagner* and Hans-Jörg Himmel Institute of Inorganic Chemistry, Ruprecht-Karls University Heidelberg, Im Neuenheimer Feld 270, 69120 Heidelberg, Germany S Supporting Information *
ABSTRACT: A new free tool for structural comparison is presented that combines existing and new features into a single software package. aRMSD incorporates the functions of establishing a pairwise correlation between the atoms of two molecular structures and the calculation of the optimal rotation matrix that minimizes the root-mean-square deviation (RMSD) between the molecules. The complexity of the Hungarian assignation problem is reduced by decomposing molecules into different subsets based on different atom or group types allowing for an efficient and robust treatment of large molecules while tolerating different substituents. Various weighting functions can be used for the calculation of RMSD values and similarity descriptors, and the utilization of coordinate uncertainties allows for the calculation of standard deviations for all calculated properties through error propagation. A new three-dimensional (3D) graphical representation that combines multiple aspects of structural information is presented which is useful in the analysis of structural similarity and diversity. The capabilities of aRMSD are demonstrated by selected examples that show how the program can be utilized in the analysis of structural changes and in the correlation of structure and activity in molecules. The source code of the program can be downloaded at https://github.com/armsd/aRMSD.
■
INTRODUCTION The ongoing process of digitization and the introduction of new computer-aided technologies and automation processes over the last 20 years caused the establishment and dynamic development of the field of chemoinformatics.1,2 One of its key aspects, the determination of diversity and similarity between molecular entities and the classification of compounds based on their structural, physical, chemical, and biological properties, is today of crucial importance in a variety of fields covering all areas of the life sciences. Especially research and development in pharmaceutical and medicinal chemistry are textbook examples of the interplay between chemical and information technology and emphasizes their role in the solution of synthetic and analytic challenges.3 Based on the quantitative structure−activity relationship (QSAR), a common method to predict the affinity and selectivity of potential drugs to binding sites in proteins is to virtually screen molecules that possess a structure similar to a known active reference molecule. In addition to pure structural similarity, different properties such as electron density distributions or the electrostatic potential of molecules may be taken into account as well which lead to the development of various descriptors and scoring values to compare molecular shapes.4,5 As a result, drug design is nowadays usually a combined approach of experiment and molecular modeling the latter being used to increase the © 2017 American Chemical Society
efficiency of the process and to help in the selection of potential candidates.6,7 By using chemical libraries consisting of thousands of different molecules, combinatorial chemistry approaches can be used to synthesize the designated molecules in parallel which are subsequently tested for biological activity by means of high-throughput screening. This protocol dramatically reduces the time required for synthesis and analysis and increases the overall effectiveness in drug discovery and other industrial processes (e.g., enzyme catalysis) and its importance is expected to further increase in the near future. In addition to the generation of libraries and databases, molecular similarities are also often used to investigate and quantify the dynamic behavior of macromolecules (e.g., trajectories from molecular dynamics (MD) simulations) in which the change of a given similarity descriptor over time can be directly linked to the flexibility or rigidity of a single molecular fragment or an entire protein domain. Furthermore, when experimental structures are compared with the results from in silico modeling, similarities can be used to estimate the quality of the given model or computational level of theory and help in the development and evaluation of new and existing simulation methods.8,9 The information obtained from Received: September 1, 2016 Published: February 13, 2017 428
DOI: 10.1021/acs.jcim.6b00516 J. Chem. Inf. Model. 2017, 57, 428−438
Article
Journal of Chemical Information and Modeling
As of now, there is a very limited number of tools available that combine the solutions to the pairwise correlation and leastsquares superposition problems. Recently, Allen and Rizzo reported the implementation of a routine to determine the optimal correspondence RMSD in the DOCK package to improve the calculation of docking success rates of symmetric molecules.31−33 A related approach combining graph theory and structural overlap was developed by Helmich and Sierka for similarity detection in global optimization problems.34 In this contribution, we present a new program package called aRMSD (automatic RMSD) which is a useful and efficient tool that matches two structures, calculates their optimal weighted RMSD, and conveniently analyzes similarity and diversity between molecules. Some of its key features are • User friendliness: Starting from the loading of files, the user is guided through the processes of matching and aligning in an interactive fashion through clearly structured menus. Whenever possible, tasks and calculation steps are carried out automatically to minimize user input and all results are directly printed on screen alongside plain and concise instructions, warnings and error messages. • File support: The program is equipped with an internal parser for some of the most common file types (xyz, mol/sdf, mol2, and out files from the computational chemistry programs Gaussian35 and Orca 36) and supports additional file formats through an interface to openbabel.37 Structural data from X-ray experiments (cif and SHELX formats) is supported as well and passed on to a menu that allows for an interactive manipulation and expansion of the structure which is afterward converted from fractional to Cartesian coordinates. Furthermore, aRMSD has its own file type (xyzs) which contains information about atomic coordinates including standard deviations and additional plotting information. • Flexibility: Atoms that belong to the same group in the periodic table are recognized and supported in the matching procedure and different weighting functions for the RMSD calculation are readily available. In addition, almost all calculation variables and default settings are loaded from an external file and can be individually adjusted by the user. • Calculation results: Besides the total RMSD of the superposed molecules and other descriptors, the contribution of individual atom types and user-defined substructures are automatically calculated. The contributions of bond distances, angles, and dihedral angles can be estimated with the help of z-matrices and all internal coordinates as well as selected individual properties can be directly and conveniently compared. • Error propagation: Standard deviations of nuclear coordinates (e.g., in crystallographic data) are automatically taken into account in all calculations which allows for an in-depth comparison between experimental structures and provides information about the statistic significance of the results as well as the uncertainties of all calculated properties (e.g., RMSD and R2 values, bond angles, etc.). • Graphical output: Publication-quality figures of different graphical representations can be generated at any stage of the program by the VTK38 software system.39 In addition to the classic superposition, a new representation of the
structural comparisons is not only useful in systematic investigations but also provides insight into the statistic significance of results. In light of the manifold applications, various scoring coefficients and descriptors for similarity recognition and have been developed over the years,10 yet the most common and straightforward method to determine the similarity between two molecular structures, defined by their coordinate vectors P and Q, is to calculate the root-mean-square deviation (RMSD, eq 1) between the interatomic distances of the equivalent atomic coordinates.11 RMSD(P, Q) =
1 N
∑ ∥Q i − Pi∥2 i
(1)
This approach is suited for macromolecules and small compounds alike, and the RMSD has been shown to be a powerful descriptor that is frequently used in structural biology.12−16 The general approach is to minimize its value by a combination of rigid body translations to a common origin and the rotation of one molecule onto the other, and several solutions to this problem have been developed over the last 50 years.17 The optimal RMSD value for two given molecules can for instance be obtained by the Kabsch algorithm.18,19 Starting from the covariance matrix A calculated from the two coordinate vectors (eq 2), the rotation matrix U which aligns the structures and minimizes their RMSD in the process can be determined by singular value decomposition (SVD, see the Supporting Information for further details about the relation of the SVD and the determination of the optimal rotation matrix) of A, as shown in eqs 3 and 4. A = PTQ
(2)
VSWT = SVD(A)
(3)
⎛1 0 0⎞ ⎜ ⎟ U = W⎜ 0 1 0 ⎟VT ⎜ ⎟ ⎝0 0 d ⎠
d = sign(det(WVT)) (4)
Alongside the Kabsch algorithm, other methods for the solution of the superposition problem were described, most notable the quaternion approaches, based on the work of Kearsley and others, which are computationally efficient and numerically stable.20−26 However, all methods have in common that a correct one-to-one correspondence between the atoms in the two molecules must exist because the optimization of the RMSD assumes equivalent atoms in both structures. If this prerequisite is not given and the coordinate sequences are not identical, they need to be established prior to the RMSD calculation. For very small systems the matching problem may be solved by hand, however even for small to medium-sized molecules this process is already unfeasible which makes it desirable to determine the correct pairwise correlation automatically. Although this problem can be solved by different methods, the Hungarian algorithm is the most widely used approach.27−29 This algorithm calculates an initial cost matrix from various properties such as interatomic distances or relative orientations of the atoms in the molecules and provides the optimal solution of the assignment problem through an iterative depletion of the matrix.30 The application of the Hungarian algorithm to larger molecules can however become computationally demanding if the data is unstructured and the full coordinate sets are used due to its O(n4) scaling. 429
DOI: 10.1021/acs.jcim.6b00516 J. Chem. Inf. Model. 2017, 57, 428−438
Article
Journal of Chemical Information and Modeling
Figure 1. Protocol of the simplified program flow of aRMSD in its present form, structured into three main aspects.
■
Step IFile Handling and Structural Consistency. The first step carried out by the program is to read-in structural data (atom types and nuclear coordinates) for two molecules from different files. File types directly supported by aRMSD are xyz, mol2, and out files from Gaussian and Orca,48 and data from additional formats can be used via an interface to openbabel. If the investigation involves X-ray structures, cifs49 can be parsed to utilize the standard deviations of the nuclear coordinates in all subsequent calculations. In these cases, the molecules may require an additional treatment (e.g., generation of symmetry equivalent atoms, removal of solvent molecules), and the fractional coordinates can be expanded and modified in an interactive fashion by the user based on an arbitrary combination of symmetry operations.50 After the full molecules have been constructed, the coordinates and their standard deviations are transformed to Cartesian coordinates. Furthermore, aRMSD reads and writes information from/to its own open file type (xyzs) which is structured into different sections and contains information about atom types, nuclear coordinates, standard deviations, bonds, and various plotting-related aspects such as colors, radii, and the camera viewpoint. By
superposed structures was developed which combines absolute deviations between individual atoms and their relative contributions to the total RMSD. By exporting data to xyzs files, it is possible to quickly pass the results to an external 3D renderer (e.g., Blender40).
COMPUTATIONAL DETAILS
aRMSD is written in Python41,42 and uses numpy43 and Cython44,45 for improved computational efficiency. Error propagation is handled by the uncertanties46 package; VTK is used for 3D rendering, and matplotlib,47 for 2D plots. All data generated by aRMSD can be exported to formatted text files which can be interfaced to third-party programs for plotting and further modification. aRMSD was designed as a Python module but can also be compiled into a standalone executable from its source code which is available at https://github.com/ armsd/aRMSD. The basic functional principle of the program is shown in Figure 1 and consists of three main aspects which will be discussed below in more detail. 430
DOI: 10.1021/acs.jcim.6b00516 J. Chem. Inf. Model. 2017, 57, 428−438
Article
Journal of Chemical Information and Modeling
Figure 2. Schematic representation of the transformation of two molecules (taken from the second example) from their initial orientation to standard orientation by translating the center of mass and rotating the eigenvectors (ωi) of the diagonalized moment of inertia tensor onto the Cartesian x, y, and z axes.
Figure 3. Decomposition of the molecules into N subsets based on atom types or groups in the periodic table in order to reduce the computational effort of the matching procedure.
storing data in xyzs files, it is possible to directly reuse crystallographic data without carrying out structural modifications and to transfer51 the data to other 3D rendering platforms such as Blender if a respective interface is provided. In order for the Kabsch algorithm to work properly, the molecular structures must be identical with respect to the number of atoms and a correct pairwise correlation between the atoms of the molecules must be established. Especially if experimental and computational data is compared, it is often the case that this prerequisite is not met as for instance hydrogen atoms may have been removed for the X-ray structure depiction which itself functions as the source for structural comparisons. Additional problems can arise due to disorder and indistinguishable atoms that formally occupy52 the same position and a correct pairwise correlation between the atoms is usually only given if the experimental structures are used as computational input. To facilitate the use of data from different file formats and to address these potential problems, aRMSD performs checks of the internal data structures and tries to automatically establish consistency between the molecules aiming to remove the need to manipulate the files by hand prior to their use. The first step is to test whether the molecules have the same number of atoms and if two or more atoms occupy the same position. If different atomic numbers are found, the program tries to establish consistency through subsequent removal of hydrogen atoms bound to carbon, to other carbon group elements and finally by removing all hydrogen atoms in the entire molecules. Occupation problems have to be solved interactively for each occurring case and the user can select the atoms from a list of possible choices. During this process the atom types are determined as they are used
during the establishment of a pairwise correlation in the matching process. If different elements are found which belong at least to the same group in the periodic table (e.g., matching of chloro and bromo substituents), the molecules are considered suitable for the execution of the Hungarian algorithm and the coordinate unit is checked by analyzing the number of bonds and converted to angstroms if needed. After the establishment of structural consistency, the data is passed to the next routine. Step IIAlignment and Matching. The prerequisite of the Kabsch algorithm is that the coordinate vectors P and Q have an identical pairwise sequence which is in most cases not given. As the matching processes relies on absolute distances (and/or relative orientations) between the potential atom pairs of the two molecules, a sufficiently good initial alignment between the structures is required. The starting point for this procedure are the molecules in their standard orientation in which the center of mass is shifted to the origin of the Cartesian coordinate system and the rotation axes ωi obtained from the diagonalized moment of inertia tensor are rotated onto the principal coordinate axes, as shown in Figure 2.53 The alignment in standard orientation can be improved by a series of arbitrary symmetry operations (rotations around the x, y, and z axes, reflections at the xy, xz, and yz planes and an inversion at the Cartesian origin) that are applied to the first of the two molecules (“Model”). After a sufficient similarity54 has been established, the atoms in the two molecules are structured, according to Figure 3, into subsets based on different atom types (e.g., H, C, O, ...) or groups (e.g., halogens, chalcogens, ...). For each subset i, a cost matrix Ci (eq 5) is calculated using a function f, that accounts for the absolute and/or relative 431
DOI: 10.1021/acs.jcim.6b00516 J. Chem. Inf. Model. 2017, 57, 428−438
Article
Journal of Chemical Information and Modeling
Figure 4. Run time of the Hungarian algorithm with (black) and without (gray) subset decomposition for molecules of different sizes with a 1:2:3:6:8 element ratio. The inset shows that the additional decomposition is the determining factor for smaller molecules up to ≈200 atoms.
calculation type is the unweighted RMSD in which all atoms are treated equally; however, in many cases it is useful to assign different weights do individual atoms. Atomic properties such as masses, nuclear charges, or interpolated structure factors56 for a given source energy can be used to acknowledge the fact that the positions of lighter atoms are experimentally not as accurately determined as they are for heavier atoms. After the determination of the optimal rotation matrix and the alignment of P and Q, various descriptors of the superposition quality are calculated, one of the most basic being the weighted (root-)mean-square deviation (R)MSD. The contributions of different structure parts (e.g., a userdefined set of atoms or individual atom types) are calculated as well which allows to investigate local contributions of fragments or domains to the superposition of the molecules in detail. Descriptors such as R2 values, cosine similarities, and GARD57 values can be calculated for the (sub)structures and a decomposition of the RMSD into internal coordinates (distances, bond angles, and dihedral angles) can be achieved by the generation and comparison of z-matrices for both molecules. Furthermore, all bond distances, average distances per bond type, bond angles, and dihedral angles are compared automatically (see first example) upon request and the parameters with the highest deviations are displayed to the user. An interpolation between the structures can be carried out as well which may serve as input in minimum energy crossing point (MECP) and transition state optimizations. Traditionally the direct superposition of the global structure or a selected fragment is shown to indicate the overall quality of the results. In addition to the classic representation, approaches that combine multiple sources of information into one single graphical representation have been proposed in the literature to further increase its meaningfulness.58 To include both absolute RMSD values associated with each atom pair as well as their relative contribution to the total RMSD of the superposition, we developed a new graphical representation (see second and
distances between all atoms of the subset and the optimal solution to each cost matrix is obtained by the Hungarian algorithm. ⎛ f (P1, Q ) ⋯ f (P1, Q ) ⎞ 1 N 1N ⎟ ⎜ 11 ⎟ ⋮ ⋱ ⋮ Ci = ⎜ ⎟⎟ ⎜⎜ ⎝ fN1 (PN , Q 1) ⋯ fNN (PN , Q N )⎠
(5)
In some cases, errors in the matching results can occur due to high differences between the full structures or a few internal coordinates (e.g., dihedral angles) which cannot be resolved by any symmetry operation and have to be manually adjusted prior to the calculation of the optimal rotation matrix. Technical Performance of the Matching Algorithm. The performance of the matching procedure was tested on a notebook with an Intel i7-3520 M 2.90 GHz processor and 8 GB RAM using molecules consisting of up to 6000 atoms belonging to five different atom types in a 1:2:3:6:8 ratio. The scaling behavior is given in Figure 4 which shows that up to ca. 200 atoms no speed-up is gained due to subset decomposition. This behavior is expected because of the additional decomposition step which outweighs the reduction of the cost matrix size. Since the longer run-time is negligible for small molecules, subset decomposition should always be used because it reduces the number of possible mismatches between the two molecular structures and ensures that only atoms of the same element or group ar correlated. Around 200 atoms the benefits from smaller cost matrices surpass the overhead of the decomposition and the efficiency is increased by a factor of ≈2.3 and allows for the matching of 15 000 atoms in less than 10 s.55 Step IIISuperposition and Postprocessing. After the establishment of the correct one-to-one correspondence between the atoms, the Kabsch algorithm is used to minimize the RMSD between the two molecules. The most common 432
DOI: 10.1021/acs.jcim.6b00516 J. Chem. Inf. Model. 2017, 57, 428−438
Article
Journal of Chemical Information and Modeling
Figure 5. Results of the direct comparison of bond distances, bond angles, and dihedral angles between the X-ray structure of [H2B(hpp)]2 and the energy minimum obtained at the SCS-MP2/def2-TZVPP level of theory.
(e.g., DFT functionals, basis sets) and the evaluation of existing ones. Furthermore, the screening of different methods is a useful and necessary step to determine an appropriate model which can afterward be applied to a whole set of related chemical problems. As part of our studies regarding the dehydrocoupling of boranes and diboranes,59,60 we attempted to benchmark different computational methods in order to determine the best approach for the calculation of structural characteristics and thermodynamical properties. Since data for the latter is scarce, the starting point was to evaluate the similarity between our computations and the experimental structure of different boranes, including the base-stabilized diborane [H2B(hpp)]2 whose Lewis structure can be seen in Figure 5.61 A complete molecule of the experimental structure62 (space group P21/c) was created with the internal cif parser of aRMSD using a 1 − x, −y, 1 − z symmetry operation and was afterward converted to Cartesian coordinates. Using the symmetrized X-ray structure as input, several DFT (density functional theory) and wave function based methods were tested for their capability to reproduce the experiment. The quality of all computational methods, summarized in Table 1, is very good which can be attributed to the almost perfect C2h symmetry of [H2B(hpp)]2 in the solid state which restricts the number of degrees of freedom and reduces the amount of possible conformers. The systematic investigation allows to differentiate between the influence of the basis set and the correlation method, and it can be seen that by going from a smaller def2-SV(P) to a larger def2-TZVP basis set the agreement between theory and experiment is improved, additional polarization functions (def2-TZVPP) however have almost no beneficial effect. Within the range of DFT methods it can be seen that hybrid functionals are slightly better than the BP86 reference. The best results however are obtained by the two MP2 approaches which perform equally within the statistic significance of the results (slightly better results are obtained for the SCS-MP2 approach if uncertainties of the experimental structure are neglected) and the superposition with the SCS-
third example) which uses the structure of the second of the two molecules (“Reference”) shown in a modified ball-and-stick model. The sizes of the sphere radii at the nuclear coordinates are proportional to the relative contribution of the respective atom pair to the total RMSD. The usage of different sphere sizes is particularly useful in cases of very small deviations between the two structures. The absolute RMSD value of each pair is given by the sphere color on a red-yellow-green (RYG) scale with a default range from ≥0.7 Å to 0.0 Å. In addition the RMSD values of each atom pair, information about differences in internal coordinates are given as well. Corresponding bonds in the molecules that are shorter or longer than a user-defined threshold (default: |0.02| Å) are shown with an intersection in green (shorter bond in the model) or red (shorter bond in the reference). In this representation (named “aRMSD representation” in the following) up to four atoms can be selected to return the RMSD value of the selected atom pair or to compare a distance, angle or dihedral angle in the two structures. The combination of absolute and relative contributions as well as internal coordinates combines multiple structural aspects in a single graphical representation which is useful to quickly unravel differences between molecules, as will be shown in the following examples.
■
RESULTS AND DISCUSSION To demonstrate the capabilities and performance of the present program, various experimentally and theoretically determined molecular structures were matched and compared by the procedures described above. Three examples were chosen to show possible applications of aRMSD in different research areas. A more comprehensive introduction and additional background information on the presented chemistry can be found in the cited literature. Influence of Computational Models on the Results of Structural Optimizations. In addition to the accuracy of relative energies, the quality of optimized structures are of key importance in the development of new computational methods 433
DOI: 10.1021/acs.jcim.6b00516 J. Chem. Inf. Model. 2017, 57, 428−438
Article
Journal of Chemical Information and Modeling
as a model for C−H and C−C bond activation processes which are thought to be important intermediates in catalysis but are difficult to study in the very labile alkane−metal σ-complexes.63 In base-stabilized diboranes two principal bonding interactions, shown left in Figure 6, are conceivable that involve different degrees of B−H and B−B contributions which were shown to depend on the nature of the transition metal fragment coordinated to the B2H2 moiety.64 Besides vibrational and NMR spectroscopy, structural characterization and analysis is of key importance to rationalize the bonding situation. We reported the synthesis and characterization of the complexes [M(CO)4{HB(hpp)}2] (M = Cr, Mo, or W) which are obtained by photolysis of the hexacarbonyls in the presence of the diborane. To compare the crystallographic data, the carbonyl complexes were parsed directly from their cif files65 into the program and the coordinates of M = Mo and W (space group C2/m) were expanded by an x, − y, z operation in order to generate the symmetry equivalent atoms along the mirror plane perpendicular to the B−B bond axis. Figure 6 shows the results of the superposition between the [HB(hpp)]2 fragment of the complex and the free diborane as reference. Neglecting differences in the orientation of CH2 groups of the bicyclic guanidines, the major changes are found in the B 2 H 2 substructure as expected. Upon coordination, the B−H bond lengths increase by an average of 0.08(3) Å while the B−B bond is shortened by 0.031(4) Å indicating a primarily transfer of electron density from the B−H bonds to the metal center which is supported by a strong redshift of the B−H vibration and quantum-chemical bonding analysis. Since the investigation of bonding interactions relies to a high degree on computational approaches, we again compared the modeling results (idealized C2v symmetric structures) with the experiment in order to ensure the quality of the structures used in our theoretical work and the results are summarized in Table 2. The overall agreement between the computational model and the experimental reference is very good in all cases as can be seen from the low RMSD and R2 and cosine similarity values close to unity.
Table 1. Quality of the Optimized Structure of the BaseStabilized Diborane [H2B(hpp)]2 with Respect to the Level of Theory (Correlation Method and Basis Set) Determined by the RMSD from the Experimental Structure Obtained by X-ray Diffractiona level of theory
def2-SV(P)
def2-TZVP
def2-TZVPP
BP86 B3LYP BHLYP PBE0 TPSSH MP2 SCS-MP2
0.102(4) 0.101(4)
0.094(4) 0.094(4)
0.094(4) 0.093(4) 0.092(4) 0.090(4) 0.090(4) 0.085(4) 0.085(4)
a
All values are calculated without the contributions of the hydrogen atoms from the CH2 groups of the bicyclic guanidines. The number of significant digits and the standard deviations originate from the quality of the experimental structure.
MP2/def2-TZVPP optimization is discussed in more detail. The contribution of the four hydrogen atoms in the molecules to the RMSD is 69.1% which is not unexpected since their positions can not be accurately determined by X-ray diffraction experiments due to the absence of core electrons. Consequently if all hydrogen atoms are excluded, the RMSD is reduced to 0.0625(6) Å. A systematic analysis between the experimental structure and the optimization is given in Figure 5. An almost ideal correlation between theory and experiment is found for all different internal coordinates, and the highest deviations for bond distances, bond angles, and dihedral angles are ≈0.1 Å, ≈6.6°, and ≈9.3°, respectively, and involve the hydrogen atoms of the BH2 groups. The results not only reveal that the structures obtained by the wave function methods are better than their DFT counterparts but also show the limits of the comparison due to the low precision of the hydrogen atom positions in the diffraction experiment. Bonding in Transition Metal Complexes of Base-Stabilized Diboranes. We have been interested in the coordination chemistry of base-stabilized diboranes due to the fact that these molecules are isoelectronic to alkanes and can therefore be used
Figure 6. Extreme cases of bonding interactions between base-stabilized diboranes (left) which can be classified as a combination of B−H−M and B−B−M centered and aRMSD representation (right) the experimental structures of [Mo(CO)4{HB(hpp)}2] (diborane fragment) and [HB(hpp)]2 which shows the changes within the B2H2 moiety upon complexation. 434
DOI: 10.1021/acs.jcim.6b00516 J. Chem. Inf. Model. 2017, 57, 428−438
Article
Journal of Chemical Information and Modeling
the atom type decompositions in Table 2. The figure also shows how individual properties are analyzed in the aRMSD representation and the atoms C-6 and C-7 are selected which directly outputs the bond distances in both molecules and their difference to the user on screen. Valence Tautomerism of Cu-GFA Complexes. Valence tautomers are compounds with two stable isomers that differ only by their electron distribution and can be interconverted by an intramolecular electron transfer induced by an external stimulus such as temperature or pressure and are of high interest due to the distinct differences of their physical properties. Recently a CuCl2 complex with a redox-active pyridine66 based guanidino-functionalized aromatic67,68 ligand (GFA, complex 2 in Figure 8, middle) that shows valence tautomerism in acetone solution was reported by us in which the fully reversible redox process involves an electron transfer from the neutral GFA ligand to the CuII centers. This reaction causes the 2-fold oxidation of the redox-active ligand to GFA2+ and the reduction of the copper centers to CuI.69 Interestingly the related complexes of the benzene70,71 or naphthalene72 based GFAs (Figure 8 left and right) do not exhibit valence tautomerism. The question why only complex 2 is able to switch between the electronic states is difficult to answer because of the similar electronic structures and ligand properties of the GFAs. To explain the observations, we analyzed the structural changes that occur upon an intramolecular electron transfer which interconverts between the high spin (triplet) and low spin (singlet) states for all three complexes.73 Optimizations at the BP86/def2-SVP level of theory without symmetry constraints were found to yield good results with respect to the experimental structures and the different spin states were compared after the exclusion of all hydrogen atoms. Figure 9 shows their superpositions in aRMSD representation, and it can directly be seen that the two structures of the pyridine system are almost identical while there are distinct structural differences between the two spin states of the other complexes and the respective RMSD values are 0.224 Å (1, 81.6% on Cl), 0.084 Å (2, 53.7% on Cl), and 0.515 Å (3, 81.7% on Cl). The fact that the structural
Table 2. Results for the Superposition of the Molecular Structures of [M(CO)4{HB(hpp)}2] (M = Cr, Mo, or W) Obtained by the X-ray Diffraction Experiment with the Results of an Optimization at the BP86/def2-SV(P) Level of Theory property
M = Cr
M = Mo
M=W
RMSD/Å wRMSDa/Å R2 cosine similarity substructure
0.2187(7) 0.1709(6) 0.99591(3) 0.99671(3) RMSD/Å
0.107(1) 0.0880(9) 0.99905(2) 0.99953(2) RMSD/Å
0.097(1) 0.074(1) 0.99922(2) 0.99967(2) RMSD/Å
[M(CO)5] 0.0898(7) 0.1327(6) [{HB(hpp)}2] 0.2504(9) 0.096(2) atom type RMSD/Å (%) RMSD/Å (%)
0.131(1) 0.081(2) RMSD/Å (%)
B (no. 2) C (no. 18) H (no. 2) N (no. 6) O (no. 4) M (no. 1)
0.039(2) (01.33) 0.2848(6) (70.23) 0.10(2) (09.40) 0.0818(9) (05.79) 0.112(1) (10.88) 0.0524(6) (02.38)
0.048(1) (03.16) 0.1025(5) (14.14) 0.16(1) (33.70) 0.0509(5) (03.49) 0.1701(9) (38.95) 0.0693(2) (06.46)
0.046(3) (04.25) 0.1006(9) (20.09) 0.06(3) (07.40) 0.0491(9) (04.78) 0.154(2) (47.06) 0.09091(3) (16.41)
a Structure factors for Mo Kα radiation were used in the calculation of the weighted RMSD.
The higher RMSD values for M = Cr can be directly related to the orientation of a single CH2 group (59.1% of the total RMSD) in the crystal structure which adopts a conformation that breaks the C2v symmetry of the molecule. Structure factors were calculated using atomic scattering factors with an X-ray source energy of 17479.34 eV (Mo Kα radiation) and used to show the effect of weights on the results. Due to the increased significance of the metal centers which already had an excellent structural correlation (low RMSD values in the decomposition per atom type in Table 2), the wRMSD is significantly lower compared to the unweighted case. The graphical results for the chromium complex are shown on the left side in Figure 7 while the aRMSD representation of the tungsten analogue is shown on the right side. The latter confirms the even distribution of the RMSD over the entire molecule which is also reflected by
Figure 7. From left to right: direct superposition (wire-frame model, left) and aRMSD representation (middle) of the structures of [Cr(CO)4{HB(hpp)}2] obtained by X-ray diffraction and optimization at the BP86/def2-SV(P) level of theory. The high contribution of a single CH2 group to the RMSD can be directly seen in both depictions. The results of the analogue complex [W(CO)4{HB(hpp)}2] with two selected atoms is shown on the right side. Hydrogen atoms bound to carbon were excluded in all calculations. 435
DOI: 10.1021/acs.jcim.6b00516 J. Chem. Inf. Model. 2017, 57, 428−438
Article
Journal of Chemical Information and Modeling
Figure 8. Lewis structures of the CuCl2 complexes (high spin states) of the GFA ligands with different aromatic cores. Only the pyridine complex 2 shows valence tautomerism in solution.
Figure 9. aRMSD representations (triplet vs singlet states) of the CuCl2 complexes of three different GFA ligands, aromatic core from left to right: benzene (RMSD = 0.224 Å), pyridine (RMSD = 0.084 Å), and naphthalene (RMSD = 0.515 Å). Hydrogen atoms were excluded in all calculations.
a more efficient investigation of macromolecules, the reduction of matching errors through better cost matrices, and the implementation of new features which will be added to aRMSD.
deviations in complex 3 are located on one side of the molecule can be attributed to the unsymmetrical molecules (see the Supporting Information). In addition, the RMSD in the pyridine complex is evenly distributed over the molecule, the contributions of the chlorido ligands are significantly lower compared to the other compounds and the number of bond lengths that differ by more than |0.02| Å is the smallest one in the set. The barrier for an electron transfer event is expected to decrease with increasing similarity between the structures of the two electronic states according to Marcus theory and the computational results predict the smallest barrier for the pyridine based ligand which is in line with the experimental observations. All three examples show how structural superpositions can help to explain experimental observations and indicate the broad range of the possible applications in various research fields.
■
ASSOCIATED CONTENT
* Supporting Information S
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jcim.6b00516. Details about the Kabsch algorithm, additional information about technical aspects of the program, and computational details (PDF)
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected].
■
ORCID
CONCLUSIONS Although the foundations of superpositions between molecular structures were formulated half a century ago, they continue to play an important role in a variety of different established and innovative research fields. Here we present a new, free userfriendly and flexible tool for the comparison of molecular structures whose central functions are a combination of the optimal atomic matching by the Hungarian algorithm and the calculation of the optimal rotation matrix based on the solution of Kabsch. In addition some unique features are implemented such as the direct support of crystallographic data and the inclusion of uncertainties through error propagation. Using the rendering capabilities of VTK, it is possible to generate graphical representations of the calculation results and a new depiction is presented which combines different aspects of structural information. The capabilities of aRMSD were demonstrated by three examples that show the importance of detailed structural analysis in the solution of chemical problems. Current development focuses on an improved renderer to allow
Arne Wagner: 0000-0002-4232-7015 Notes
The authors declare no competing financial interest.
■ ■
ACKNOWLEDGMENTS The authors thank the Deutsche Forschungsgemeinschaf t (DFG) for continuous financial support. REFERENCES
(1) Engel, T. Basic Overview of Chemoinformatics. J. Chem. Inf. Model. 2006, 46, 2267−2277. (2) Agrafiotis, D. K.; Bandyopadhyay, D.; Wegner, J. K.; van Vlijmen, H. Recent Advances in Chemoinformatics. J. Chem. Inf. Model. 2007, 47, 1279−1293. (3) Willett, P. Similarity Methods in Chemoinformatics. Annu. Rev. Inf. Sci. Technol. 2009, 43, 1−117. (4) Maldonado, A. G.; Doucet, J.; Petitjean, M.; Fan, B.-T. Molecular Similarity and Diversity in Chemoinformatics: From Theory to Applications. Mol. Diversity 2006, 10, 39−79. 436
DOI: 10.1021/acs.jcim.6b00516 J. Chem. Inf. Model. 2017, 57, 428−438
Article
Journal of Chemical Information and Modeling (5) Nikolova, N.; Jaworska, J. Approaches to Measure Chemical Similarity. QSAR Comb. Sci. 2003, 22, 1006−1026. (6) Sliwoski, G.; Kothiwale, S.; Meiler, J.; Lowe, E. W., Jr. Computational Methods in Drug Discovery. Pharmacol. Rev. 2014, 66, 334−395. (7) Kapetanovic, I. Computer-Aided Drug Discovery and Development (CADDD): In Silico-Chemico-Biological Approach. Chem.-Biol. Interact. 2008, 171, 165−176. (8) Grimme, S. A General Quantum Mechanically Derived Force Field (QMDFF) for Molecules and Condensed Phase Simulations. J. Chem. Theory Comput. 2014, 10, 4497−4514. (9) Grimme, S.; Brandenburg, J. G.; Bannwarth, C.; Hansen, A. Consistent Structures and Interactions by Density Functional Theory with Small Atomic Orbital Basis Sets. J. Chem. Phys. 2015, 143, 054107. (10) Todeschini, R.; Consonni, V.; Xiang, H.; Holliday, J.; Buscema, M.; Willett, P. Similarity Coefficients for Binary Chemoinformatics Data: Overview and Extended Comparison Using Simulated and Real Data Sets. J. Chem. Inf. Model. 2012, 52, 2884−2901. (11) McLachlan, A. D. A Mathematical Procedure for Superimposing Atomic Coordinates of Proteins. Acta Crystallogr., Sect. A: Cryst. Phys., Diffr., Theor. Gen. Crystallogr. 1972, A28, 656−657. (12) Cohen, F. E.; Sternberg, M. J. On the Prediction of Protein Structure: The Significance of the Root-Mean-Square Deviation. J. Mol. Biol. 1980, 138, 321−333. (13) Brüschweiler, R. Efficient RMSD Measures for the Comparison of two Molecular Ensembles. Proteins: Struct., Funct., Genet. 2003, 50, 26−34. (14) Mechelke, M.; Habeck, M. Robust Probabilistic Superposition and Comparison of Protein Structures. BMC Bioinf. 2010, 11, 363. (15) Gapsys, V.; de Groot, B. L. Optimal Superpositioning of Flexible Molecule Ensembles. Biophys. J. 2013, 104, 196−207. (16) Hung, L.-H.; Samudrala, R. Accelerated Protein Structure Comparison using TM-score-GPU. Bioinformatics 2012, 28, 2191− 2192. (17) Diamond, R. A Mathematical Model-Building Procedure for Proteins. Acta Crystallogr. 1966, 21, 253−266. (18) Kabsch, W. A Solution for the Best Rotation to Relate two Sets of Vectors. Acta Crystallogr., Sect. A: Cryst. Phys., Diffr., Theor. Gen. Crystallogr. 1976, A32, 922−923. (19) Kabsch, W. A Discussion of the Solution for the Best Rotation to Relate two Sets of Vectors. Acta Crystallogr., Sect. A: Cryst. Phys., Diffr., Theor. Gen. Crystallogr. 1978, A34, 827−828. (20) Kearsley, S. K. On the Orthogonal Transformation Used for Structural Comparisons. Acta Crystallogr., Sect. A: Found. Crystallogr. 1989, A45, 208−210. (21) Diamond, R. A Note on the Rotational Superposition Problem. Acta Crystallogr., Sect. A: Found. Crystallogr. 1988, A44, 211−216. (22) Horn, B. K. P. Closed-form Solution of Absolute Orientation Using Unit Quaternions. J. Opt. Soc. Am. A 1987, 4, 629−642. (23) Theobald, D. L. Rapid Calculation of RMSDs Using a Quaternion-Based Characteristic Polynomial. Acta Crystallogr., Sect. A: Found. Crystallogr. 2005, A61, 478−480. (24) Kneller, G. R. Superposition of Molecular Structures Using Quaternions. Mol. Simul. 1991, 7, 113−119. (25) Coutsias, E. A.; Seok, C.; Dill, K. A. Using Quaternions to Calculate RMSD. J. Comput. Chem. 2004, 25, 1849−1857. (26) Liu, P.; Agrafiotis, D. K.; Theobald, D. L. Fast Determination of the Optimal Rotational Matrix for Macromolecular Superpositions. J. Comput. Chem. 2009, 31, 1561−1563. (27) Kuhn, H. W. The Hungarian Method for the Assignment Problem. Nav. Res. Logist. Q. 1955, 2, 83−97. (28) Kuhn, H. W. Variants of the Hungarian Method for Assignment Problems. Nav. Res. Logist. Q. 1956, 3, 253−258. (29) Munkres, J. Algorithms for the Assignment and Transportation Problems. J. Soc. Ind. Appl. Math. 1957, 5, 32−38. (30) The Hungarian algorithm is one of the basic algorithms which can be used to solve the assignation problem, i.e., to find an optimal assignment between two sets of variables. The correspondence
between the two variable sets is represented by a cost matrix whose matrix element ij correspond to the penalty between the elements i and j. Through a series of repeated subtractions (and additions) of the smallest elements from all rows and columns of the entire matrix the optimal assignment is obtained if the number of lines required to cover all zero elements is equal to the number of the entries in the variable sets. (31) Lang, P. T.; Brozell, S. R.; Mukherjee, S.; Pettersen, E. F.; Meng, E. C.; Thomas, V.; Rizzo, R. C.; Case, D. A.; James, T. L.; Kuntz, I. D. DOCK 6: Combining Techniques to Model RNA-Small Molecule Complexes. RNA 2009, 15, 1219−1230. (32) Brozell, S. R.; Mukherjee, S.; Balius, T. E.; Roe, D. R.; Case, D. A.; Rizzo, R. C. Evaluation of DOCK 6 as a Pose Generation and Database Enrichment Tool. J. Comput.-Aided Mol. Des. 2012, 26, 749− 73. (33) Allen, W. J.; Rizzo, R. C. Implementation of the Hungarian Algorithm to Account for Ligand Symmetry and Similarity in Structure-Based Design. J. Chem. Inf. Model. 2014, 54, 518−529. (34) Helmich, B.; Sierka, M. Similarity Recognition of Molecular Structures by Optimal Atomic Matching and Rotational Superposition. J. Comput. Chem. 2012, 33, 134−140. (35) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö .; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, Revision B.01; Gaussian Inc.: Wallingford CT, 2010. (36) Neese, F. The ORCA Program System. WIREs Comput. Mol. Sci. 2012, 2, 73−78. (37) O’Boyle, N. M.; Banck, M.; James, C. A.; Morley, C.; Vandermeersch, T.; Hutchison, G. R. Open Babel: An Open Chemical Toolbox. J. Cheminf. 2011, 3, 33. (38) Schroeder, W.; Martin, K.; Lorensen, B. The Visualization Toolkit, 4th ed.; Kitware, 2006. (39) The 3D rendering capabilities are optimized to produce highquality images which makes the visualization of large molecules with a few 1000 atoms difficult on notebooks. (40) Blender Online Community. Blender−a 3D Modelling and Rendering Package; Blender Foundation, Blender Institute: Amsterdam, 2016. (41) Python Software Foundation. Python Language Reference, version 2.7/3.5. Available at http://www.python.org (accessed Dec 28, 2016). (42) Oliphant, T. E. Python for Scientific Computing. Comput. Sci. Eng. 2007, 9, 10−20. (43) van der Walt, S.; Colbert, S. C.; Varoquaux, G. The NumPy Array: A Structure for Efficient Numerical Computation. Comput. Sci. Eng. 2011, 13, 22−30. (44) Behnel, S.; Bradshaw, R.; Citro, C.; Dalcin, L.; Seljebotn, D.; Smith, K. Cython: The Best of Both Worlds. Comput. Sci. Eng. 2011, 13, 31−39. (45) Bradshaw, R.; Behnel, S.; Seljebotn, D. S.; Ewing, G. The Cython Compiler. http://cython.org (accessed Dec 28, 2016). (46) Lebigot, E. O. Uncertainties: a Python Package for Calculations with Uncertainties. http://pythonhosted.org/uncertainties/ (accessed Dec 28, 2016). (47) Hunter, J. D. Matplotlib: A 2D Graphics Environment. Comput. Sci. Eng. 2007, 9, 90−95. 437
DOI: 10.1021/acs.jcim.6b00516 J. Chem. Inf. Model. 2017, 57, 428−438
Article
Journal of Chemical Information and Modeling (48) If information about the Cartesian gradient exist in the files, it can be used as “standard deviation” in subsequent calculations by request of the user. (49) Hall, S. R.; Allen, F. H.; Brown, I. D. The Crystallographic Information File (CIF): A new Standard Archive File for Crystallography. Acta Crystallogr., Sect. A: Found. Crystallogr. 1991, A47, 655−685. (50) Depending on the structural data, some insight concerning the crystal symmetry and the operations is required to construct a complete molecule from the initial fragments may be required from the user. (51) A script for automatic rendering of data from xyzs files with Blender is part of the aRMSD release. (52) In case of indistinguishable atoms (e.g., C and N atoms) usually both are written to file with the same nuclear coordinates. (53) All graphical representations of molecular structures have been created with the internal VTK renderer of aRMSD. (54) At this stage of the algorithm, “sufficient similarity” corresponds to a qualitatively correct orientation for both molecules in which the relevant structural aspects (e.g., rings, metal-fragments, side chains) are properly aligned. The solution provided by the standard orientation approach has to be accepted by the user prior to the matching procedure and or can be modified by symmetry operations. (55) In the same time roughly 10000 atoms can be matched by the Hungarian algorithm without subset decomposition. (56) Henke, B. L.; Gullikson, E. M.; Davis, J. C. X-Ray Interactions: Photoabsorption, Scattering, Transmission, and Reflection at E = 50− 30,000 eV, Z = 1−92. At. Data Nucl. Data Tables 1993, 54, 181−342. (57) Baber, J. C.; Thompson, D. C.; Cross, J. B.; Humblet, C. GARD: A Generally Applicable Replacement for RMSD. J. Chem. Inf. Model. 2009, 49, 1889−1900. (58) Scott, W. R. P.; Straus, S. K. Determining and Visualizing Flexibility in Protein Structures. Proteins: Struct., Funct., Genet. 2015, 83, 820−826. (59) Schulenberg, N.; Ciobanu, O.; Kaifer, E.; Wadepohl, H.; Himmel, H.-J. The Doubly Base-Stabilized Diborane(4) [HB(μhpp)] 2 (hpp = 1,3,4,6,7,8-hexahydro-2H-pyrimido[1,2-a]pyrimidinate): Synthesis by Catalytic Dehydrogenation and Reactions with S8 and Disulfides. Eur. J. Inorg. Chem. 2010, 2010, 5201−5210. (60) Wagner, A.; Litters, S.; Elias, J.; Kaifer, E.; Himmel, H.-J. Chemistry of Guanidinate-Stabilised Diboranes: Transition Metal Catalysed Dehydrocoupling and Hydride Abstraction. Chem. - Eur. J. 2014, 20, 12514−12527. (61) Ciobanu, O.; Allouti, F.; Roquette, P.; Leingang, S.; Enders, M.; Wadepohl, H.; Himmel, H.-J. Thermal and Catalytic Dehydrogenation of the Guanidine-Borane Adducts H3B·hppH (hppH = 1,3,4,6,7,8hexahydro-2H-pyrimido[1,2-a]pyrimidine) and H3B·N(H)C(NMe2)2 A Combined Experimental and Quantum Chemical Study. Eur. J. Inorg. Chem. 2008, 2008, 5482−5493. (62) CCDC reference number: 690038. (63) Wagner, A.; Kaifer, E.; Himmel, H.-J. Diborane(4)-Metal Bonding: Between Hydrogen Bridges and Frustrated Oxidative Addition. Chem. Commun. 2012, 48, 5277−5279. (64) Wagner, A.; Kaifer, E.; Himmel, H.-J. Bonding in DiboraneMetal Complexes: A Quantum Chemical and Experimental Study of Complexes Featuring Early and Late Transition Metals. Chem. - Eur. J. 2013, 19, 7395−7409. (65) CCDC reference numbers: 865135 (M = Cr), 865133 (M = Mo), and 865132 (M = W). (66) Stang, S.; Lebkücher, A.; Walter, P.; Kaifer, E.; Himmel, H.-J. Synthesis of Redox-Active Guanidine Ligands with Pyridine and pBenzoquinone Backbones. Eur. J. Inorg. Chem. 2012, 2012, 4833− 4845. (67) Himmel, H.-J. Guanidinyl-Functionalized Aromatic Compounds (GFAs) − Charge and Spin Density Studies as Starting Points for the Development of a New Class of Redox-active Ligands. Z. Anorg. Allg. Chem. 2013, 639, 1940−1952.
(68) Himmel, H.-J. Redox-Active Guanidines and GuanidinateSubstituted Diboranes. Top. Heterocycl. Chem. 2015, DOI: 10.1007/ 7081_2015_168. (69) Wiesner, S.; Wagner, A.; Kaifer, E.; Himmel, H.-J. A Valence Tautomeric Dinuclear Copper Tetrakisguanidine Complex. Chem. Eur. J. 2016, 22, 10438−10445. (70) Peters, A.; Kaifer, E.; Himmel, H.-J. 1,2,4,5-Tetrakis(tetramethylguanidino)benzene: Synthesis and Properties of a New Molecular Electron Donor. Eur. J. Org. Chem. 2008, 2008, 5907−5914. (71) Peters, A.; Trumm, C.; Reinmuth, M.; Emeljanenko, D.; Kaifer, E.; Himmel, H.-J. On the Chemistry of the Strong Organic Electron Donor 1,2,4,5-Tetrakis(tetramethylguanidino)benzene: Electron Transfer in Donor-Acceptor Couples and Binuclear Late Transition Metal Complexes. Eur. J. Inorg. Chem. 2009, 2009, 3791−3800. (72) Vitske, V.; König, C.; Kaifer, E.; Hübner, O.; Himmel, H.-J. Syntheses of the First Coordination Compounds of the New Strong Molecular Electron Donor and Double Proton Sponge 1,4,5,8Tetrakis(tetramethylguanidino)naphthalene. Eur. J. Inorg. Chem. 2010, 2010, 115−126. (73) Wiesner, S.; Wagner, A.; Kaifer, E.; Himmel, H.-J. The Control of the Electronic Structure of Dinuclear Copper Complexes of RedoxActive Tetrakisguanidine Ligands by the Environment. Dalton Trans. 2016, 45, 15828−15839.
438
DOI: 10.1021/acs.jcim.6b00516 J. Chem. Inf. Model. 2017, 57, 428−438
