{'fraction_non_alphanumeric': 0.03789341001299878, 'fraction_numerical': 0.014259266553747981, 'mean_word_length': 4.702605570530099, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Revisiting van der Waals Radii: From Comprehensive Structural Analysis to Knowledge‐Based Classification of Interatomic Contacts†

Abstract Weak noncovalent interactions are responsible for structure and properties of almost all supramolecular systems, such as nucleic acids, enzymes, and pharmaceutical crystals. However, the analysis of their significance and structural role is not straightforward and commonly requires model studies. Herein, we describe an efficient and universal approach for the analysis of noncovalent interactions and determination of van der Waals radii using the line‐of‐sight (LoS) concept. The LoS allows to unambiguously identify and classify the “direct” interatomic contacts in complex molecular systems. This approach not only provides an improved theoretical base to molecular “sizes” but also enables the quantitative analysis of specificity, anisotropy, and steric effects of intermolecular interactions.


Introduction
Van der Waals (vdW) radii (R vdW ) are one of the cornerstones of contemporary chemistry. The visual perception and basic analysis of chemical structures and architectures are strongly related to the concept of vdW radii, which provides the basic definition of the atomic "sizes". Even though this concept has little physical grounding, it is easily interpretable and therefore widely used in the theoretical chemistry, particularly in design of force fields (in the Lennard-Jones parameter form), [1] dispersion corrections, [2] and COSMO-RS-type solvate models. [3] These models operate -directly or indirectly -by the intermolecular interatomic distances, which are regarded as corresponding to the minimum energy of atom-atom interactions, and are usually defined as a function of the sum of vdW radii. However, it has been noted [4,5] that the sum of R vdW for available vdW sets including the most widely used one tabulated by Bondi [6] consistently underestimate the position of the energy minimum by as much as 0.3-0.4 Å. An illustrative example is a C 2 Cl 6 crystal with the shortest r(Cl ··· Cl) = 3.65 Å that exceeds significantly the 2R vdW B°ndi (Cl) = 3.50 Å. The latter can only be achieved at a very high pressure of ca. 1.2 GPa. [7] Moreover, about 2.0 % (> 4000) of unique organic molecular crystals in the Cambridge Structural Database (CSD) [8] do not contain intermolecular contacts shorter than the sum of R vdW B°ndi (see the first section of Supporting Information for details). The question arises: why are all bonds in these crystals longer than the distance, which is usually considered as the most probable, in other words, closest to the energy minimum of the respective interatomic potential?
This inconsistency is rooted in the indirect nature of the approaches used to obtain the vdW parameters from experimental structural datasets. Conventionally, R vdW are derived from the position of the vdW peak in the distributions of contact distances between nonbonded atoms (Figure 1a,b). [6,9,10] Such distributions represent a superposition of a gaussian curve due to the vdW peak, and a rapidly growing function due to the randomly distributed contacts (Figure 1a). Because the position of the vdW peak maximum (D max ) is often hidden within the dataset for other contacts (Figure 1b), the sum of vdW radii is normally taken as the half-height of the vdW peak distance, D half instead of the D max ( Figure 1a). However, D max has a clearer physical meaning as corresponding to the most probable distance for the considered interaction. [4] On the contrary, D half -based vdW radii are qualitative, as the interpretation of their sum is not straightforward. It should be mentioned here that the relation between statistical distributions and energy surfaces is not straightforward and cannot be quantitatively determined without introducing arbitrary and unwarranted assumptions. [11] However, low-energy regions of molec- [a] I. Y. Chernyshov, Prof. Dr. E. A. Pidko ular potential energy surfaces can be recognized and mapped from distributions of observed structures, thus D max can be considered an approximation of the "minimum energy position". [4,11] Although alternative knowledge-based [12] approaches to derive R vdW were proposed, similar to the original Bondi's strategy they fail to extract the most probable interatomic distances from the structural data. In addition to statistical approaches, several ab initio computational ways to derive atomic sizes and R vdW have been proposed. [13,14] However, the real systems are much more complex than the computational models. As a result, such computational approaches provide only qualitative radii values and their use for the analysis of real extended systems therefore necessitates further investigations to ensure an adequate connection with the experimental results.
In this work, we introduce a new method to define the D max and, accordingly, to compute D max -based vdW radii (R max ) for almost any intermolecular interaction directly from diverse structural datasets. The key idea is to eliminate the background from randomly distributed contacts using the Line-of-Sight [15] concept. It considers a pair of atoms as interacting only when they "see" each other (Figure 1d) because no other atom in the structure intrudes between them or, in other words, shields them from one other. This concept allows one to redetermine vdW parameters from the bulk structural data and provide an opportunity to find atom-type specific D max -based vdW radii for elements in different chemical environments with clear physical meaning. It should be noted that intermolecular interactions cannot be reduced to atom-atom interactions [16] and energy decomposition schemes should rather be applied for a detailed analysis. [17] Nevertheless, our approach is in fact a statistical way to localize the most probable geometries of intermolecular interactions, which is independent of the interpretation of the results in terms of atom-atom interactions. We argue that R max can be used to improve the accuracy of computational approaches directly or indirectly involving consideration of the size of molecules. It will be useful for the analysis of intermolecular interactions in combination with quantum chemistry methods.

Results and Discussion
We propose to implement the Line-of-Sight [15] concept (LoS) to eliminate the background from randomly distributed contacts, which is equivalent to classification of contacts into corresponding and noncorresponding to interatomic interactions. Thus, we start with understanding the physical meaning of this procedure. LoS concept considers the atoms as interacting only when they "see" each other (Figure 1d) because no other atom in the structure intrudes between them or, in other words, shields them from one another (for the exact definition see the SI). We identify three main types of contact shielding, namely, the covalent, intra-and intermolecular shielding depending on the type of shielding atom (Figure 1d), with only the covalent shielding having the characteristic distances comparable to those of the D max . A representative example of the contributions from the LoS and shielded contacts for the complete dataset of N ··· O contacts is illustrated in Figure 1c.
To qualitatively validate and illustrate this concept, we carried out a detailed conformational and bonding analysis on a model H 3 N ··· CH 3 F system ( Figure 2) by means of density functional theory (DFT) calculations (B3LYP/6-311 + + G(d,p)). The ab initio electron density-based noncovalent interaction surface (NCIS) method [18] was employed to detect and estimate the intermolecular interactions. Basic geometric considerations within the LoS model imply the N ··· C interaction for a configuration with ffN ··· CÀ F of 180°and ffN ··· CÀ H of~70°. The decrease of ffN ··· CÀ F would result in the gradual shielding of the N ··· C contact by the H atoms and ultimately vanishing the N ··· C interaction at ffN ··· CÀ F~155-160°and ffN ··· CÀ H~45-50°, depending on the assumed R vdW . The NCIS analysis supports these predictions ( Figure 2). The NCIS has a pronounced maximum between the N and C atoms in linear H 3 N ··· CH 3 F. Upon bending the geometry, the NCIS redistributes towards the N ··· H area evidencing strengthening of the respective N ··· H contact with a concomitant weakening of the initial N ··· C interaction ( Figure 2). These data show that the shielding does not fully eliminate the original two-atom contact but rather decreases its contribution to the overall intermolecular interaction. On contrary, the LoS contacts are dominated by their respective diatomic contributions that makes them perfectly suitable for the statistical analysis of the intermolecular contacts. Note that strictly speaking the LoS approach cannot be used as a criterion of significance of intermolecular interactions. The presented results and analysis are not sufficient to unambiguously evaluate whether only LoS contacts give rise to significant intermolecular interactions.
The LoS concept was next used to determine the vdW radii for the main group elements from the structural data in the Cambridge Structural Database (CSD). [8] We processed 224 001 unique CSD entries, from which more than 40 000 000 intermolecular LoS contacts were obtained and analyzed. An iterative procedure was employed, in which at the 0 th iteration, the LoS was defined with Bondi's R vdW . [6] The LoS contact distance distributions were used to determine D max values for all possible A ··· B atom pair contacts, while the D max values were used to compute the vdW radii, R vdW , for all elements in the dataset by means of least-squares minimization of the function [Eq. (1)]: where the outer summation is carried out over all selected A ··· B contacts, and w(A ··· B) is a weight defined as the root square of the number of LoS contacts shorter than D max . This weight function was chosen to eliminate bias caused by the determination of D max values of distributions based on a small number of contacts from the vdW radii. However, the use of the degenerate weight function (w(A ··· B) = 1) gives almost the same R vdW values with the largest deviation of 0.02 Å. Because the contact shielding depends on vdW radii, the R vdW values obtained at the first step were used to recalculate the contact shielding parameters and determine the D max values for all atom pairs. This procedure was repeated iteratively until convergence of 0.001 Å was reached for all R vdW . For some atom pairs, it was not possible to accurately locate D max due to the imperfect shape of the gaussian peak in the LoS distribution. This can be attributed to one of the following three scenarios, namely, (1) a small total number of contacts, (2) the analyzed contact type consists of several contact subtypes characterized by different D max values, resulting in a significant broadening of the vdW peak, and (3) significant dependence of the D max on the contact geometry, which also widens the vdW peak. Addressing the second scenario by splitting the contact type into subtypes reduces the number of contacts and therefore often gives rise to the first scenario. As a result, we were able to obtain D max for type-specific contacts with C (C sp2 and C sp3 ) and O (O sp2 and O sp3 ) atoms only. On the contrary, for B, P, and As we were not able to accurately obtain D max for any contact type. It should be noted that for N, S, and Se, which have several common atom types, we were able to identify one atom type, namely, the N sp3 , S sp3 and Se sp3 , as other types of these elements exhibit a significant dependence on contact geometry (the second scenario).
Applying the R vdW search procedure to the selected contact types has revealed that for some contacts D max significantly deviates from the sum of the atomic type-specific R vdW (Table S1). Such deviations are observed either when the contact corresponds to a potentially specific interaction (e. g. hydrogen or halogen bond) or the interacting atoms are possibly sterically hindered. Such contacts were not used to determine R vdW as our final goal was to obtain R vdW corresponding to weak nonspecific interactions that were not affected by secondary effects. The contacts used for R vdW determination are marked with ' + ' sign in the second column of Table S1. The R vdW values determined using this procedure (R max ) are listed in Table 1 and compared with the respective values from the most popular vdW radii sets (R half ).
The results in Table 1 reveals that the LoS model yields vdW radii consistently exceeding the values obtained by D half -based approaches by as much as 0.10-0.20 Å resulting in the increase of the sum of R vdW by 0.2-0.4 Å. This is attributed to the transition of the analysis from the rather arbitrary D half parameters to more specific and rigorously defined D max values. We argue that R max are more physically sound as they directly reflect the most probable contact distance, whereas R half appear to correspond to strongly shortened interactions, which are usually denoted as "specific" in the chemical literature.
Interestingly, the vdW radii estimates from quantum chemical calculations for isolated atoms (R QM ) by Rahm et al., [14] defined as the average distance from the nucleus to a point where the electron density falls to 0.001 a.u., also exceed R half (Table 1). The comparison with the data obtained in this work shows that these theoretical R QM values cohere well with R max for C, N, O, and F, whereas they quite uniformly exceed R max for S, Se, Cl, Br and I by 0.15-0.20 Å. This implies that the electron density parameters for weak interactions at the most probable distance depend only weakly on the interacting elements.
Furthermore, our data analysis allows identifying the contact-angle dependencies in vdW radii. A representative example is the interhalogen interactions CÀ Hal 1 ··· Hal 2 À C, Hal=Cl, Br, I (X-bonds, XBs), which are usually classified as types I and II. [19] Type I XBs are characterized by ffCÀ Hal 1 ··· Hal 2 � ffCÀ Hal 2 ··· Hal 1 and are usually nonspecific and weak, whereas type II XBs are characterized by ffCÀ Hal 1 ··· Hal 2 � 90°and ffCÀ Hal 2 ··· Hal 1 � 180°and are usually strong and shorter than type I XBs due to the σ-hole interaction. [ It should be noted, that R max are in fact a condensed representation of the analyzed contacts and therefore can be instrumental for the more detailed analysis of noncovalent interactions and statistical definition of qualitative structural features of chemical systems. For example, specific interactions and sterically hindered contacts can be automatically identified by D max substantially deviating from the sum of the respective vdW radii. Let us consider two representative examples of CH ··· O and C ··· O contacts, for which the D max are, respectively, shorter or longer than the sum of R max .
CH ··· O contacts show D max values shorter by 0.14 Å than the sum of R max . This deviation is much larger than that (< 0.02 Å) detected for the related CH ··· N and CH ··· F interactions (Table S1). The CH ··· O contacts are shorter by 0.18 Å than CH ··· N and very similar to CH ··· F (Figure 3a), although an opposite trend can be seen in the C sp2 ··· O/N/F contacts (Figure 3b). Such deviation of D max from the sum of R max implies that CH ··· O contacts are highly specific and their structureforming role is more significant than that of the other CH ··· X contacts, including CH ··· F, despite the similar electrostatic nature of these contact types. [22] The C CR4 ··· O contacts formed by quaternary carbon exceed by 0.30 Å the corresponding R max sum, which, in turn, cohere well with the D max values for the C ··· O contacts formed by primary (C Me ) and tertiary (C HR3 ) carbon atoms (Figure 4). Such a shift of the C CR4 ··· O vdW peak is in line with the expected significant steric repulsions between the neighbors of the interacting carbon atom and the oxygen atom and indicates significant steric hinderance of C CR4 atoms. This means that the concept of vdW radii should be applied with a special care when dealing with the atoms in confined tetrahedral or octahedral environments (e. g. Si, Bi, most of d, f-elements, and R max Bondi [6] R.&T. [9] Alv. [10] H other elements considered by Alvarez in Ref. [10]), which would show an increased effective size due to the unaccounted steric effects. The corresponding D max values will merely indicate the position of the first maximum of the radial distribution function and make little sense in terms of noncovalent interactions. Only the atoms in an "open" environment such as trigonal, square planar or square pyramidal configurations can be used to determine R vdW from contact distance distributions without shielding effects to be crucial. Even for the extended dataset considered herein these conditions were fully satisfied only for the 10 elements, which data are summarized in Table 1. These examples illustrate that despite all advantages, the element-defined R max radii still fail to generally define the most probable distances of intermolecular interactions. The preferred and more accurate approach is to directly utilize D max values determined for the given atom pair contact rather than the sum of vdW radii. In this way, the structural analysis will automatically account for the specificity of the analyzed interactions, their anisotropy, and the impact of steric effects. We have developed an algorithm and implemented it in a script (deposited at the GitHub [23] ) that allows searching the CSD for intermolecular contacts with certain geometry (ConQuest output) with their subsequent classification as line-of-sight or shielded contacts ( Figure 5). The produced datasets can readily be used for the determination of the D max using any standard table processing or data analysis software. Such D max values obtained for specific contacts with a certain geometry can be  [21] is used instead of histograms for the purpose of clarity.  used to create high-throughput geometry-based descriptors of intermolecular interactions, applicable e. g. to high-throughput screening of heterogeneous catalysts. [24] Although the described approach is statistical, it can be used for the analysis of specific systems. In this case, one needs to compare the length of an intermolecular contact with D max of the distribution of the corresponding contacts. Detailization of the environment of interacting atoms and interaction geometry eliminates the bias caused by differences in the nature of atoms and anisotropy of the interaction. As a result, one can localize the most probable geometry of the specific condensed system using experimental structural data. The same outcome can be achieved by analyzing potential energy surface by means of quantum chemical modeling, however, its application to condensed systems is qualitative and not straightforward. [25] As a representative example of the approach we provide the analysis of the short F ··· F contact in crystalline pentafluorobenzoic acid (PFBA, CSD refcode PFBZAC01). [26] It is 2.63 Å long, and both CÀ F ··· F are equal to 155°. This is much shorter than the sum of the element-specific R max (3.10 Å) suggesting that such a contact should be forced and potentially repulsive in nature. However, this conclusion does not hold when atomtype specific radii are utilized for the analysis. The interatomic contact cannot anymore be regarded as an ultrashort one when considering D max for linear F ··· F contacts (2.90 Å), or, even more specifically, the aromatic C 3 F 3 fragments (2.78 Å). (Figure 6). Unfortunately, further refinement of the fluorine-containing fragments is currently not possible due to the small amount of available structural data. However, given that the fluorine environment in PFBA is more electronegative than in more common C 6 F 5 groups, we can assume that the most probable r (F ··· F) in PFBA is even smaller than 2.78 Å. This implies that such "ultrashort" contact is the intrinsic characteristics of the molecule itself and the role of the packing features for their occurrence is only minor if any.

Conclusions
In summary, a new approach for automated analysis and classification of noncovalent interactions has been developed.
The key feature of the reported methodology is the filtration of the direct interatomic contacts from the other background structural features based on the LoS concept. The validity of this approach has been confirmed by density functional theory calculations. The LoS concept was used to recompute R vdW for light elements, which turned out to be substantially larger than those currently accepted by the chemistry community. Our analysis reveals that the underestimations of the "atomic sizes" noted earlier for most tabulated vdW radii (R half ) is solely related to the deficiencies of the datasets used for their determination. The R max introduced here are free from statistical bias and are based on clear physical grounds. We anticipate the high practical utility of R max , and particularly, their atom type-specific variations for improving of molecular sizes in different chemical methods and approaches, and for the analysis of intermolecular interactions on a wide variety of systems including the experimental structural databases and the results of theoretical calculations on molecular and condensed systems relevant to various fields of chemistry and material sciences. The current LoS approach implemented now in the relevant software enables the quantitative analysis of specificity, anisotropy and steric effects of intermolecular interactions while benchmarking databases as well as studying specific systems.

Computational Methods
Version 5.39 with 4 updates (up to August 2018) of CSD [8] was used for selection of organic crystals containing H(D), B, C, N, O, F, Si, P, S, Cl, As, Se, Br and I atoms. Disordered, erroneous, polymeric, pressurized, powder structures and experiments with R-factor > 0.075 were removed from consideration. 224 001 selected crystals were used for search of unique intermolecular contacts A ··· B with distance D(A ··· B) up to 7.0 Å. CÀ H, NÀ H and OÀ H bond lengths were normalized to CCDC/ConQuest defaults: CÀ H: 1.089 Å, NÀ H: 1.015 Å, OÀ H: 0.993 Å. Acetylenic C sp À H bond lengths were normalized to 1.06 Å (neutron diffraction: ACETYL05, RALDEN01, XEHLEB, ZULDEP01), and SÀ H bond lengths were normalized to 1.34 Å (microwave data: H 2 S, [27] neutron diffraction: NALCYS02). For each contact the following information were collected: (1) contact geometry, (2) chemical nature of contact atoms, (3) shielding atom, and (4) contact shielding value. A total of c.a. 640 000 000 contacts were found, more than 40 000 000 of which was LoS depending on used vdW radii (41 346 551 for the final version from Table 1). These data were used further to build and analyze the distribution of distances of various contacts and is available for download. [28] Obtained line-of-sight A ··· B contacts were used to plot histograms of contact distance distributions for all possible A, B atom type pairs (Table S1). Primary analysis of these distributions combined with chemical common sense allowed to select atom types and, therefore, contacts used for the van der Waals radii determination. The main requirements were: (1) there are enough contacts to determine D max so that line-of-sight contacts distribution contains a line-of-sight peak that is similar to the Gaussian function; (2) contact corresponds to non-specific interaction; (3) influence of steric effects on A ··· B distance can be excluded; (4) addition of the contact to the list of contacts used for van der Waals radii Figure 6. Distance distribution for linear (CÀ F ··· F > 150°) intermolecular F ··· F contacts between two electron-withdrawing (EWG) fragments. From left to right: red dashed line denotes the shortest F ··· F intermolecular distance in PFBZAC01; orange line stands for D max of linear F ··· F contacts between two aromatic fragments; green and purple lines stand for D max of linear and all F ··· F contacts, respectively. determination does not change any radius significantly (> 0.05 Å). This selection is subjective to a certain degree, however, it's well suited for our final goal which was to obtain R vdW corresponding to weak nonspecific interactions that are not affected by any effects. Selected contact types were used for R max determination and are marked with ' + ' sign in the second column of Table S1. Using the data of Table S1 together with the data on intermolecular contacts [28] allows one to reproduce this work, or to choose another set of contact types to determine another version of R max .
Quantum chemical calculations were performed with the Gaussian 16 rev. B.01 [29] program at the B3LYP/6-311 + + G(d,p) level of theory. The NCIS analysis was performed with the Multiwfn package. [30] Geometry of the H 3 N ··· CH 3 F complex was obtained by merging optimized NH 3 and CH 3 F molecules into the staggered C 3v structure with D(N ··· C) = 3.0 Å. The N ··· CÀ F angle was changed so that the H 3 N ··· CH 3 F complex retains the N ··· (H)CÀ F symmetry plane. NCI surface (0.5 isosurface value) dependence on N ··· CÀ F angle remains unchanged with D(N ··· C) increasing.