Toward First-Principles Design of Organic Nonlinear Optical Materials

Aug 1, 2018 - NanoScience Technology Center, University of Central Florida , Orlando , Florida 32826 , United States. § Department of Chemistry, Univ...
0 downloads 0 Views 3MB Size
Article pubs.acs.org/crystal

Cite This: Cryst. Growth Des. XXXX, XXX, XXX−XXX

Toward First-Principles Design of Organic Nonlinear Optical Materials: Crystal Structure Prediction and Halogen Bonding Impact on Hyperpolarizabilities of 2‑Iodo-3-hydroxypyridine Irina D. Yushina,† Artëm E. Masunov,*,†,‡,§,∥,⊥ Diana Lopez,‡,§ Alexander A. Dyakov,† and Ekaterina V. Bartashevich† †

South Ural State University, Lenin pr. 76, Chelyabinsk 454080, Russia NanoScience Technology Center, University of Central Florida, Orlando, Florida 32826, United States § Department of Chemistry, University of Central Florida, Orlando, Florida 32816, United States ∥ Photochemistry Center RAS, Federal Research Center Crystallography and Photonics Russian Academy of Science, Ul. Novatorov 7a, Moscow 119421, Russia ⊥ National Research Nuclear University MEPhI (Moscow Engineering Physics Institute), Kashirskoye shosse 31, Moscow 115409, Russia

Downloaded via UNIV OF NORTH DAKOTA on August 1, 2018 at 12:27:29 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.



S Supporting Information *

ABSTRACT: Computational methods can potentially accelerate development of more efficient organic materials for second harmonic generation. Here, we test the method that includes the evolutionary algorithm for predicting crystal structure and prognosis of nonlinear optical properties based on the predicted structure. For this test, we selected 2-iodo-3hydroxypyridine, which exhibits second harmonic generation intensity comparable to that of urea. We performed global minimization of the lattice energy and found the experimental structure when many-body dispersion correction is added to the density functional theory values. We analyzed geometric preferences of the halogen bonding in predicted virtual polymorphs. We also found linear correlation between the lengths of the iodine−iodine halogen bonds and calculated second order susceptibilities.

I. INTRODUCTION Organic molecular1−3 and ionic4,5 nonlinear optical (NLO) materials offer a number of advantages, compared to their inorganic6,7 counterparts: ultrafast NLO response, structural versatility, and processability. Two aspects are equally important in engineering of organic NLO materials: selecting a chromophore with a large second order molecular polarizability (hyperpolarizability) and crystal packing of these chromophores in a parallel or nearly parallel noncentrosymmetric motif to maximize second order NLO response.8−10 Theoretical methods, including density functional theory (DFT), were found to be invaluable tools for computational predictions of molecular NLO properties,11−14 and hyperpolarizability in particular.15−18 Predictions of the bulk crystalline second order susceptibility χ(2) (hyperpolarizability β per unit volume) is also possible when the structure is known.19,20 Predicting the crystal packing proved to be more challenging; however, several empirical trends had been rationalized, including various supramolecular synthons21,22 and distinct patterns of hydrogen bonding.23,24 In order to compliment these trends with more predictive considerations, we recently proposed25 to employ global minimization of the lattice energy for the periodical system. A similar approach is © XXXX American Chemical Society

being used in organic crystal structure predictions (CSP), relevant to pharmacophore materials.26,27 An evolutionary algorithm, implemented in USPEX code28 can direct the search toward lower energy areas of potential surface landscape and focus a further optimization path by preserving and combining energetically favorable structural features. USPEX was already established as an efficient method for inorganic crystal structure predictions.29−31 Apart from well-known hydrogen bonds, some other important noncovalent interactions may determine both packing and NLO properties of organic crystals. These interactions include halogen, chalcogen, and pnictogen bonds, which have their own geometric preferences. These preferences were theoretically studied in the model molecular complexes.32,33 In particular, the halogen−halogen (Hal···Hal) interactions can play a key role in crystal structure formation, and affect the physical properties of crystals.34−38 Taking these interactions into account complicates the modeling, structure prediction, and crystal engineering. Yet, it may enhance the Received: April 8, 2018 Revised: July 8, 2018

A

DOI: 10.1021/acs.cgd.8b00529 Cryst. Growth Des. XXXX, XXX, XXX−XXX

Crystal Growth & Design

Article

Figure 1. (a) The halogen and hydrogen bonds in F33. (b) ELF distribution for the halogen bond I(1)···I(2); the plane is C−I(1)−I(2). (c) The criterion ρ(r)min vs φ(r)min: the electron density (black) and the electrostatic potential (red) along the halogen bond I(1)···I(2); R, interatomic distance (a.u.).

identical halogens tends to be nonpolar or pure van der Waals (vdw). For the quantitative description, one may use the features of the electron density and electrostatic potential anisotropy50 around a covalently bound halogen atom. Applying certain electronic criteria allows one to judge if a particular halogen atom in a given crystal structure plays the role of electronacceptor and to confirm that it is indeed involved in the halogen bond51 according to IUPAC definition.52 The following characteristics may be used: electron localization function (ELF),53 one-electron potential (OEP),54 Laplacian of electron density,55 superposition of gradient fields of the electron density and electrostatic potential,56 and atomic dipole polarization.57 These functions make it possible to identify the halogen bonds in complicated cases of multiple interactions58 and to distinguish between halogen bonds and weak vdw interactions in borderline cases, when the geometric criteria may not be sufficient. Anisotropy of halogen bonds presents a challenge for their modeling when empirical force fields are used. While some authors introduced vdw repulsion anisotropy in order to predict crystal structures,59 others argue that anisotropy of electrostatic interactions is sufficient and propose adding offcenter charges or atomic-centered dipole and quadrupole moments.60 The latter approach was implemented in an AMOEBA force field, which we use in this work. van der Waals parameter fitting for Cl atoms was published recently.61 These parameters were chosen to reproduce experimental densities and heats of vaporization for chloro-substituted methanes. To fit vdw AMOEBA parameters for I atoms, we used ab initio potential surfaces in model molecular dimers. Full details of this work will be published elsewhere. These new vdw parameters were applied in this work to predict crystal structure and nonlinear optical properties of 2iodo-3-hydroxypyridine (Refcode YELQUC).62 For this purpose, we followed the protocol developed in the previous work.25 This protocol includes five steps: (1) customized parameter derivation for polarizable force field AMOEBA; (2) local minimizations of trial crystal structures with these parameters, combined with an evolutionary algorithm for the global minimum search; (3) filtering out duplicate virtual polymorphs produced; (4) reoptimization and final ranking based on density functional theory (DFT) with dispersion correction; and (5) prediction of the second order crystalline susceptibility in the most stable of predicted polymorphs. We already validated this protocol in the case of urea,25 a small molecule with crystal packing controlled by hydrogen bonding.

capability of NLO materials design and accuracy in prediction of crystal structures and properties. The interplay between structure and observable NLO properties was studied before,39,40 although the particular role of noncovalent interactions has been revealed only recently.41 In some cases, NLO response can be explained by the presence of noncovalent interactions of different types. For instance, in the series of 2-amino-3-nitropiridinium halides,42 chloride and bromide with coplanar crystal packing exhibit high second harmonic generation (SHG) efficiency, while corresponding nonplanar 2-amino-3-nitropiridinium iodide is an order of magnitude less efficient. Azobenzene polymers43 present another example. There, the role of the halogen bonding is found to be more significant than that of the hydrogen bonding due to higher directionality of the halogen bond. The largest SHG efficiency was reported for iodine-substituted structure, as the iodine atom is known to be highly polarizable and an effective halogen bonding donor. Statistical analysis of Hal···Hal intermolecular contacts reveals their geometric preferences in crystals formed by halogenated molecules.44−46 These preferences become especially apparent when the structures with hydrogen bonds, which may distort molecular packing, are excluded from this analysis.34 The mutual orientations of C−Hal covalent bonds are often classified into one of two groups: Type I and Type II. In Type I contacts, two C−Hal bonds are antiparallel (either nearly so or exactly by symmetry), while in Type II contacts two covalent bonds are nearly at a right angle (L-shape geometry). Two models are often used to interpret the nature of these interactions: orbital and electrostatic ones. According to the orbital model, two types of localized molecular orbitals are considered: electron donating n (nonbonding, or lone pairs) and electron accepting σ* (CHal antibonding). The maximum overlap n−σ* of the respective localized orbitals is then used to interpret different kinds of halogen interactions. Specifically, Type II interactions are formed due to a single donor−acceptor overlap, while Type I interactions are the result of double overlap, where both Hal atoms participate as electron donors and electron acceptors.47 According to the electrostatic model, the negative region of electrostatic potential around the Hal atom, formed by the lone pairs is called the nucleophilic belt, and the positive region of electrostatic potential formed on the extension of its covalent bond is called the electrophilic σ-hole.48,49 In this model, Type II interactions are observed if the σ-hole of one halogen is directed toward the nucleophilic belt of another halogen, while the nature of Type I interactions between two B

DOI: 10.1021/acs.cgd.8b00529 Cryst. Growth Des. XXXX, XXX, XXX−XXX

Crystal Growth & Design

Article

entries in the Cambridge Structural Database (CSD). The remaining monoclinic and orthorhombic SPG were combined in five additional searches with two, four, and eight molecules in the unit cell, in order to cover over 95% of CSD entries. One independent molecule in general position was used: two molecules for SPG 02 and 04; four molecules for SPG 09, 14, 19, and 33; and eight molecules for SPG 15, 61, and 62. It is important to note that this choice of SPG only affected the initial guess. No symmetry constraints were imposed during the search, so that resulting structures often belong to a different SPG. The search starts with 48 random individual structures, each relaxed to the local minimum using the sequence of TINKER modules optrigid, optimize, and xtalmin. Redetermination of SPG is performed using the f indsym utility,75 which is a part of the USPEX package. After local optimizations are completed, 30% of individuals with higher lattice energies are discarded, and remaining ones are used in the next generation. The second and subsequent generations included the seven best individuals from the previous generation, five individuals produced with rotational mutation, and 24 individuals produced with heredity operators, as described in ref 73. The remaining 12 structures were built randomly, and local optimizations were repeated. The search continues for 90 generations and produces 48 virtual polymorphs with the lowest lattice energies in each SPG. After removal of duplicate structures (as described below), the final number of predicted polymorphs was 179. The appearance of multiple duplicate structures is expected. Even the random guesses that fall in one basin in configuration space would relax to the same local minimum. An evolutionary search increases this possibility by focusing the search on the vicinity of the low local minima that were already found. Similar structures are prevented from being used in one heredity operator but are not eliminated from the rest of the search. As a result, care should be taken to not run the search for too many generations. When this happens, all 48 top candidates are duplicates of the global minimum at the AMOEBA theory level, and structural diversity required for reranking is lost. Post processing included a packing similarity search among these 48 predictions in each SPG using the Mercury program.76 This similarity search is based on a comparison of molecular clusters, containing the central molecule and its 14 nearest neighbors. If two clusters matched within 20% of intermolecular distances and angles, the respective crystal structures were considered to belong to the same packing mode. One representative from each packing mode was reoptimized at the DFT level using the FHI-aims program (revision 160328_3).77 The PBE exchange-correlation functional, “second tier” atom-centered numerical basis set and 2 × 2 × 2 k-point integration scheme were used. Additive van der Waals energy correction by Tkatchenko and Scheffler78 was used for geometry optimization (DFT-TS theory level), while more accurate many-body dispersion (MBD) energy correction79 was employed for single point energy refinement and final stability ranking (DFT-MBD theory level). The use of multiple theory levels was determined by practical considerations. DFT-MBD theory level, including nonadditive dispersion correction, was shown to be indispensable for predicting the relative energies of molecular crystals.80 As implemented in FHI-aims rev. 160328_3, it only allows for single point energy evaluations. Therefore, structure relaxations were performed at the less advanced DFT-TS theory

Here, we turn to a more complicated case that requires iodine−iodine halogen bonds modeling.63,64 We aim to validate the protocol, identify the problems found, and offer possible solutions. The crystal structure of 2-iodo-3-hydroxypyridine (YELQUC)62 is noncentrosymmetric, belongs to space group Pna21 (SPG 33), and demonstrates nonlinear optical properties (measured SHG activity is comparable to that of urea62). Molecules of 2-iodo-3-hydroxypyridine have one rotatable bond; in the crystal, the molecules form chains by somewhat distorted N···H−O hydrogen bonds. These chains are packed in antiparallel mode and linked in a zigzag motif with the typical I···I halogen bonds (Figure 1a). In this work, we generate virtual polymorphs, rank them by lattice energy, and perform statistical analysis of noncovalent iodine’s bonds in the predicted packing modes. Next, we report analytic predictions of NLO properties in the virtual polymorphs predicted to be noncentrosymmetric.

II. COMPUTATIONAL DETAILS II.A. AMOEBA Force Field Parameters. The force field parameters were generated for the 2-iodo-3-hydroxypyridine molecule using script PolType,65 installed on XSEDE.66 This script runs the Gaussian 09 program67 and TINKER postprocessing tools68 to implemena t multistep protocol of AMOEBA force field parameter generation as described in the SI for ref 69. The geometry is optimized using the Hartree− Fock (HF) method with an all-electron 6-311G(d) basis set.70,71 The bond lengths and valence angles are used to derive valence parameters of the force field. Initial atomic charges, dipoles, and quadrupoles are derived from the wave function obtained at second order Müller-Plesset perturbation theory level MP2/6-311G(d) by means of distributed multipole analysis as implemented in the GDMA program.72 These electrostatic parameters are then refined, based on MP2/6311++G(2d,2p) electrostatic potential. Atomic dipole polarizabilities and vdw parameters are assigned based on generalized atom types. Default AMOEBA values of vdw parameters remained unchanged, except for iodine atoms. Iodine vdw parameters were set by keywords “vdwpr 401 401 4.55 0.25” (to describe I···I interactions), and “vdw 410 4.35 0.50” (to describe iodine interactions with all other atoms). Here, 401 is iodine atom type, 4.35 is vdw diameter in Å, and 0.52 is the depth of vdw minimum in kcal/mol. These parameters were selected based on ab initio potential curves for iodoethylene molecule dimers. The radial potential curves were produced for both head-on and side-on orientations with C···I and I···I contacts. In addition, angular dependence for Type I and Type II I···I contacts were studied. The benchmark interaction energies were predicted using counterpoise corrected coupled clusters with single and double substitutions and perturbative triples CCSD(T) theory level, combined with cc-pVTZ basis set and effective core potentials for I. Further details of the fitting will be published elsewhere. II.B. Search for the Most Stable Crystal Structure. The global minimization of the lattice energy was performed using an evolutionary algorithm as implemented in the program USPEX,28 version 9.4.4 in molecular crystal structure prediction mode.73 Nine separate searches were performed, one for each guess of SPG: 02 (P1̅), 04 (P21), 09 (Cc), 14 (P21/c), 15 (C2/c), 19 (P212121), 33 (Pna21), 61 (Pbca), and 62 (Pnma). These SPGs are most frequently found74 in organic crystal structures and together account for over 86% of all C

DOI: 10.1021/acs.cgd.8b00529 Cryst. Growth Des. XXXX, XXX, XXX−XXX

ÅÄÅ ÅÅ 19 d = ÅÅÅÅ ÅÅ 105 ÅÅÇ

Crystal Growth & Design level. Each evolutionary search required hundreds (if not thousands) of structure relaxations. When done at the DFT-TS theory level, it runs for several hours on 2048 processors. Such computational power was not available to us. That is why it was important to select for this purpose an inexpensive computational method, such as AMOEBA, which is still accurate enough to include the correct prediction in the top 48 candidates. II.C. Statistical Analysis of Halogen Bonds. The distribution of intermolecular I···I contacts over Type I and Type II was analyzed for three data sets. Two of them were collected from the contacts in predicted polymorph structures after optimization at the AMOEBA theory level (Tinker-set), and at the DFT-TS theory level (FHIaims-set). The third data set (CSD-set) included I···I contacts from the Cambridge Structural Database (CSD v. 5.38).81 CSD-set was formed automatically using fragment C−I···I−C, where both I atoms form covalent bonds to the aromatic C atom. The following query parameters were used: I···I distance is between 3.6 and 4.6 Å, ∠C−I···I angles θ1 and θ2 are between 0° and 180° (θ1 > θ2). Additional filters were as follows: 3D coordinates determined, no disorder, no errors, not polymeric, no ions, only organics (organometallic crystals were excluded). Each of these three data sets was used to form I···I contact distribution over |θ1 − θ2| and classified in terms of Type I and Type II interactions. II.D. Predictions of Nonlinear Optical Properties. Nonlinear susceptibility prediction and electron density analysis were performed using CRYSTAL1482 at the hybrid DFT level with the PBE0 exchange-correlation functional and Grimme dispersion correction83 (DFT-D2 theory level). The basis set DZVP84 was used for I atoms. The basis sets designed by Gatti85 were used for C, H, N, and O atoms. K-point sampling was done using a Monkhorst−Pack grid of 8 × 8 × 8. Two series of calculations were performed. In the first one, predicted crystal structures were taken as optimized by FHIaims at the DFT-TS theory level. In the second series, the experimental structure YELQUC was modified as follows. While two of the lattice constants were fixed at their experimental values, the third one scanned the range from 3.7 up to 4.5 Å, and atomic positions were optimized with CRYSTAL14. Theses variations did not change space group symmetry or overall packing motif. Static second-order polarizabilities were computed via the coupled perturbed Kohn−Sham analytical approach86 at the G point of the Brillouin zone for DFT-TS optimized structures of the predicted polymorphs. The shift to second order susceptibilities was performed according to the relation: χ(2) = (β·2π)/V

∑ (χ(2)iii )2 + i

(2)ijj +

(2)jkk

44 105

13 105

∑ χ(2)iii χ i≠j

∑ (χ(2)ijj )2 +

ÉÑ1/2 ÑÑ Ñ 5 + (χ (2)ijk )2 ÑÑÑÑ ÑÑ 7 ÑÑÖ i≠j

13 105

∑ χ(2)iij χ ijk

(3)

The experimental SHG intensity is proportional to the d2 value89 and is rather sensitive to the experimental setup. For that reason, predicted SHG intensity is often reported not in absolute values but in relative terms as a comparison to a standard SHG material, such as urea:

d2 . 2 durea

Calculations of the electron density ρ(r), the electrostatic potential φ(r), and the electron localization function (ELF)53 distributions were performed using TOPOND14.90

III. RESULTS AND DISCUSSION III.A. Crystal Structure Prediction. To predict the crystal structures for 2-iodo-3-hydroxypyridine, nine searches were performed (one for each guess space group). We started with SPG 33 (Pna21) of the guess, the symmetry known to give the correct answer. Out of 48 individual structures with low lattice energy, found within 90 generations, the packing similarity search identified seven distinct packing modes. These modes are listed as MA33, MB33, ..., MG33 in Table S1 of the Supporting Information, along with their lattice energies at various theory levels. The packing similarity search identified MF33 to match the experimentally known crystal structure of YELQUC. Among all predictions in guess SPG 33, this packing mode was ranked sixth best according to AMOEBA lattice energies. Reoptimization at the DFT-TS theory level resulted in the third best ranking for MF33, and single point DFTMBD energy evaluation identified MF33 packing mode to be the best in this search. The other eight single SPG searches (Tables S2−S9) were performed in a similar manner. They found eight more packing modes (MA09, MA14, MB14, MC14, MD14, ME14, MA19, MB19) to be more stable than the MF33 one, according to AMOEBA ranking. The DFT-TS ranking identified only three such modes (MA14, MB14, MA19), and DFT-MBD ranking found none. An additional three searches with guesses made in the remaining monoclinic SPGs (Tables S10−S12 for Z = 2, 4, 8) and two searches in remaining orthorhombic SPGs (Tables S13−14 for Z = 4, 8) after DFT-MBD ranking did not produce a polymorph candidate that is more stable or even comparable (within 0.5 kcal/mol) to the M33 one. Therefore, our search, combined with DFT-MBD ranking, found the experimentally known structure to be the global minimum of lattice energy, in agreement with the experiment. Let us consider the MF33 crystal structure in greater detail. Its crystal packing is identical to the reference structure YELQUC: molecules are linked in one-dimensional chains by O−H···N hydrogen bonds, with I···I halogen bonds between these chains. The intermolecular I···I distances and key angles θ1 and θ2 (defined in Figure 1a) are insignificant (0.1 Å and 1°). However, the O−H···N hydrogen bond deviates stronger: in MF33, the angle φ = 168.0°, whereas for the reference YELQUC, this angle is φ = 130.1°. This deviation may be due

(1)

where V corresponds to the unit cell volume. To simplify the analysis, the following averaging of the tensor elements can be introduced: vectorization87 χ(2)x = χ(2)xxx + χ(2)xyy + χ(2)xzz χ(2)y = χ(2)yyy + χ(2)xxy + χ(2)yzz χ(2)z = χ(2)zzz + χ(2)xzz + χ(2)yyz χ(2)tot = (β2x + β2y + β2z)1/2

Article

(2)

and Kurtz and Perry averaging (where i, j, k = x, y, z):88,89 D

DOI: 10.1021/acs.cgd.8b00529 Cryst. Growth Des. XXXX, XXX, XXX−XXX

Crystal Growth & Design

Article

Figure 2. Parameters of noncovalent interaction (a) according to the distribution of I···I contacts over parameter |θ1 − θ2| in (b) CSD-set, (c) Tinker-set, and (d) FHIaims-set. Type I is marked in pink. Quasi-Type I is in gray, and Type II is in blue.

Figure 3. Common “zigzag” fragments of the crystal packing that is important in exhibiting NLO properties.

arrow in Figure 1b indicates the direction of electron attraction of I(2) to the region of lowered shielding of the I(1) nucleus. One dimensional plots of electron density ρ(r) and electrostatic potential φ(r) along the line of intermolecular I···I contact is shown in Figure 1c. The order of onedimensional minima on these functions can serve as a criterion to distinguish the typical halogen bond from Type I interaction.63,91,92 The minimum ρ(r)min serves as the indicator of the interatomic boundary.93 In MF33, the atom I(2) has the

to a well-known inaccuracy in the hydrogen position as reported by X-ray diffraction. Next, we examine details of the electronic structure. A twodimensional plot of the Electron Localization Function (ELF) is presented in Figure 1b. This plot illustrates the features of the halogen bonds. For the typical halogen bond I(1)···I(2), the region with high ELF values is oriented toward the σ-hole, found on the extension of the C−I(1) covalent bond. The E

DOI: 10.1021/acs.cgd.8b00529 Cryst. Growth Des. XXXX, XXX, XXX−XXX

Crystal Growth & Design

Article

Table 1. Hydrogen and Halogen Bond Characteristics in the Natural Structure and in Five Predicted Model Structures, Optimized at DFT-TS Theory Level structure (final SPG)

type of I···I interaction II, “zigzag” II, “zigzag” II, “zigzag” I, “zigzag”

YELQUC (33) MF33 (33) MG02 (09) MG62 (05)

Type Type Type Type

MI02 (01) MO04 (04)

quasi-type I “zigzag” no I···I contacts

I(1)···I(2), Å

θ1, deg

θ2, deg

|θ1 − θ2|, deg

N···H, Å

φ(OHN), deg

3.738 3.750 3.836 4.053 4.231 4.631 4.150

174.5 177.0 164.1 136.9 119.6 116.3 137.6

114.5 108.7 100.2 123.2 106.9 116.3 107.5

60.0 68.3 64.3 13.7 12.7 0 30.1

2.120 1.594 1.583 1.733

130.1 169.4 175.3 174.1

1.608 1.782

172.3 174.4

Figure 4. Intermolecular interactions in optimized structures predicted to be noncentrosymmetric.

φ(r)min minimum closer to its nucleus (Figure 1c), and it is the electron donor atom because its fraction of electrons is attracted to the I(1) nucleus strictly along the σ-hole direction. Generally, for the Type I interaction, the positions of ρ(r)min and φ(r)min minima coincide with each other. III.B. Statistical Analysis. Since our prediction produced a sizable collection of the virtual polymorphs, it would be interesting to analyze the distribution of intermolecular I···I contacts in those structures and compare their characteristics to those in naturally occurring C(sp2)−I···I−C(sp2) interactions in molecular crystals. The respective histograms are shown in Figure 2. Following Desiraju et al.,44 we indicated regions on histograms which correspond to Type I, Type II, and the intermediate “Quasi-Type I.” Distribution of geometric parameters of I···I interactions at CSD-set shows that for I atoms in molecular crystals, Type II interaction prevails over Type I. This result agrees with the earlier statistical analysis,44 cited in the review.51 Predicted structures in Tinker-set and FHIaims-set show that, respectively, there are 33% and 67% more Type I contacts than Type II ones, while in CSD-set, there are 7% more Type II contacts than Type I ones. As we

can see in the histograms (Figure 2) and Tables S15−S17, the number of I···I interactions, comparing Tinker-set and FHIaims-set, has changed qualitatively and quantitatively. Comparing Tinker-set with FHIaims-set, one can see that the second maximum of distribution shifted toward the Type I region and the overall number of I···I interactions increased. Consequently, we can conclude that either the molecular structure, presence of stronger intermolecular interactions (Hbonding), or DFT-TS optimization increases the number of Type I interactions at the expense of Quasi-Type I interactions. III.C. Second Order Susceptibility Prediction. Finally, we estimate second order NLO characteristics among predicted crystal structures and elucidate the cumulative influence of I···I noncovalent interactions and hydrogen bonds on NLO properties. Here, we consider the predicted structures that belong to five noncentrosymmetric SPGs: P1 (SPG 01), P21 (SPG 04), C2 (SPG 05), Cc (SPG 09), and Pna21 (SPG 33), which demonstrated the lowest DFT-MBD energy within their SPG. These structures show a variety of packing styles (Figure 3): different types of “zigzag” motifs based on I···I interactions (MF33, MG02, MG62, MI02) and F

DOI: 10.1021/acs.cgd.8b00529 Cryst. Growth Des. XXXX, XXX, XXX−XXX

Crystal Growth & Design

Article

MF33 demonstrates only one nonzero vectorial component χ(2)Z, which coincides with the I···I zigzag chain in the crystal. One may hypothesize that shorter I···I distances would amplify halogen bonding interaction and increase NLO response. To validate this hypothesis, we performed a series of structure optimizations with fixed a and b lattice parameters and a varied c parameter in the range from 3.7 to 4.5 Å (see section II.D for details). All molecules retained their orientation in the unit cell, indicating that the change of one lattice parameter value does not change the overall packing style of MF33. We found nearly linear dependence of calculated βtot value on c length (Figure 5c). Similar dependence was found when we varied a and b parameters (Figure 5a,b). Next, we plotted NLO properties vs I···I distances in the crystalline structures with modified lattice parameters. The dependence of I···I distances vs β (Figure 6)

no I···I interactions at all (MO04). Hydrogen bonds retain the structure-forming interaction, being almost equally short in all the cases considered (Table 1). The interaction I···I turns out to be less prevalent. For instance, in the analyzed subset, we can see packing with the absence of I···I short contacts (MO04), with the Type I contact (MG62), and with the Quasi-Type I contact (MI02). Type II contact (halogen bonds) exists in the crystal structure MF33 that corresponds to experimental crystal packing as well as in MG02 with similar geometric characteristics. The largest values of χ(2)tot and d are observed for the structure MG02 (SPG 09) with a “zigzag” motif formed by halogen bonds; this χ(2)tot value is almost identical to the one for MF33 structure (SPG 33). In addition to I···I contacts, one may consider OH···N hydrogen bonds. They link molecules in the chains along the a axis in MG02 and MF33 structures (Figure 4). In the case of MG02, the hydrogen bonded chains are parallel, while in MF33 they are antiparallel. This difference is affecting the χ(2)X component of susceptibility and results in higher overall values of χ(2)tot for MG02, although the χ(2)Z component in this case is lower than that in MF33 and the halogen bond is also longer in comparison to MF33 (Table 2). Table 2. Calculated NLO Characteristics in Five Predicted Model Structures Optimized at the DFT-TS Theory Level SPG χ(2)tot nonzero contribution d/durea

MG02

MF33

MO04

MI02

MG62

09 38.33 χ(2)Z, χ(2)X 5.37

33 34.70 χ(2)Z

04 24.05 χ(2)Y

05 0.01 χ(2)Y

5.03

3.65

01 14.59 χ(2)X, χ(2)Y, χ(2)Z 3.12

Figure 6. Relationship between βtot values and length of I···I halogen bond in the row of structures with artificially changed unit cell lattice parameters for a, b, and c crystallographic axes. Star corresponds to the optimized reference structure YELQUC.

0.36

remains linear, revealing that this descriptor is the primary reason for an NLO increase or decrease. The hydrogen bond length does not demonstrate a similar effect on NLO properties (Table S18 and Figure S2).

Polymorph MO04 has a very interesting structure without I···I interactions, but with a similar orientation of aromatic molecules. This example shows the possibility of noncentrosymmetric molecular packing even without stabilizing I···I interactions, thus halogen bonds in the case of 2-iodo-3hydroxypyridine is not the only way to direct the crystal packing. This structure demonstrates a significantly lower χ(2)tot value; however, in MI02, a quasi-type I zigzag with distant I···I interaction results in even lower values of χ(2)tot. Finally, in the considered set of structures almost planar packing in MG62 with antiparallel orientation of the molecules results in negligible values of χ(2)tot. Thus, I···I interaction of Type II plays an important role in NLO properties of 2-iodo-3hydroxypyridine. This conclusion is supported by the fact that

IV. CONCLUSIONS We used a combination of the USPEX evolutionary algorithm with AMOEBA force field local structure optimizations to predict the crystal structure of 2-iodo-3-hydroxypyridine. Noncentrosymmetric molecular packing in this structure appears to be determined by intermolecular I···I halogen bonds. The predicted crystal structure MF33 matched the experimentally observed one. It was found to be ninth best according to AMOEBA ranking, fourth best according to DFT-

Figure 5. Relationship between βtot values and artificially changed cell axis length for a, b, and c crystallographic axes. In all cases, hollow dots correspond to the optimized reference structure YELQUC. G

DOI: 10.1021/acs.cgd.8b00529 Cryst. Growth Des. XXXX, XXX, XXX−XXX

Crystal Growth & Design

Article

(2) Castet, F.; Rodriguez, V.; Pozzo, J. L.; Ducasse, L.; Plaquet, A.; Champagne, B. Design and characterization of molecular nonlinear optical switches. Acc. Chem. Res. 2013, 46 (11), 2656−2665. (3) Ratera, I.; Veciana, J. Playing with organic radicals as building blocks for functional molecular materials. Chem. Soc. Rev. 2012, 41 (1), 303−349. (4) Cui, Y. J.; Li, B.; He, H. J.; Zhou, W.; Chen, B. L.; Qian, G. D. Metal-organic frameworks as platforms for functional materials. Acc. Chem. Res. 2016, 49 (3), 483−493. (5) Liu, X.; Yang, Z.; Wang, D.; Cao, H. Molecular structures and second-order nonlinear optical properties of ionic organic crystal materials. Crystals 2016, 6 (12), 158. (6) Rondinelli, J. M.; Kioupakis, E. Predicting and Designing Optical Properties of Inorganic Materials. In Annual Review of Materials Research; Clarke, D. R., Ed.; Annual Reviews, 2015; Vol. 45, pp 491− 518. (7) Petrov, V.; Ghotbi, M.; Kokabee, O.; Esteban-Martin, A.; Noack, F.; Gaydardzhiev, A.; Nikolov, I.; Tzankov, P.; Buchvarov, I.; Miyata, K.; Majchrowski, A.; Kityk, I. V.; Rotermund, F.; Michalski, E.; Ebrahim-Zadeh, M. Femtosecond nonlinear frequency conversion based on BiB3O6. Laser Photonics Rev. 2010, 4 (1), 53−98. (8) Wampler, R. D.; Begue, N. J.; Simpson, G. J. Molecular design strategies for optimizing the nonlinear optical properties of chiral crystals. Cryst. Growth Des. 2008, 8 (8), 2589−2594. (9) Ostroverkhov, V.; Singer, K. D.; Petschek, R. G. Secondharmonic generation in nonpolar chiral materials: relationship between molecular and macroscopic properties. J. Opt. Soc. Am. B 2001, 18 (12), 1858−1865. (10) Suponitsky, K. Y.; Masunov, A. E. Supramolecular step in design of nonlinear optical materials: Effect of pi···pi stacking aggregation on hyperpolarizability. J. Chem. Phys. 2013, 139 (9), 094310. (11) Masunov, A.; Tretiak, S.; Hong, J. W.; Liu, B.; Bazan, G. C. Theoretical study of the effects of solvent environment on photophysical properties and electronic structure of paracyclophane chromophores. J. Chem. Phys. 2005, 122 (22), 224505. (12) De Boni, L.; Toro, C.; Masunov, A. E.; Hernandez, F. E. Untangling the excited states of DR1 in solution: An experimental and theoretical study. J. Phys. Chem. A 2008, 112 (17), 3886−3890. (13) Mikhailov, I. A.; Bondar, M. V.; Belfield, K. D.; Masunov, A. E. Electronic properties of a new two-photon absorbing fluorene derivative: the role of hartree-fock exchange in the density functional theory design of improved nonlinear chromophores. J. Phys. Chem. C 2009, 113 (48), 20719−20724. (14) Belfield, K. D.; Bondar, M. V.; Hernandez, F. E.; Masunov, A. E.; Mikhailov, I. A.; Morales, A. R.; Przhonska, O. V.; Yao, S. Twophoton absorption properties of new fluorene-based singlet oxygen photosensitizers. J. Phys. Chem. C 2009, 113 (11), 4706−4711. (15) Suponitsky, K. Y.; Tafur, S.; Masunov, A. E. Applicability of hybrid density functional theory methods to calculation of molecular hyperpolarizability. J. Chem. Phys. 2008, 129 (4), 044109. (16) Suponitsky, K. Y.; Masunov, A. E.; Antipin, M. Y. Computational search for nonlinear optical materials: are polarization functions important in the hyperpolarizability predictions of molecules and aggregates? Mendeleev Commun. 2009, 19 (6), 311−313. (17) Suponitsky, K. Y.; Liao, Y.; Masunov, A. E. Electronic hyperpolarizabilities for donor-acceptor molecules with long conjugated bridges: calculations versus experiment. J. Phys. Chem. A 2009, 113 (41), 10994−11001. (18) Suponitsky, K. Y.; Masunov, A. E.; Antipin, M. Y. Conformational dependence of the first molecular hyperpolarizability in the computational design of nonlinear optical materials for optical switching. Mendeleev Commun. 2008, 18 (5), 265−267. (19) Draguta, S.; Fonari, M. S.; Masunov, A. E.; Zazueta, J.; Sullivan, S.; Antipin, M. Y.; Timofeeva, T. V. New acentric materials constructed from aminopyridines and 4-nitrophenol. CrystEngComm 2013, 15 (23), 4700−4710. (20) Yushina, I. D.; Batalov, V. I.; Bartashevich, E. V.; Davydov, A. O.; Zelenovskiy, P. S.; Masunov, A. E. Raman spectroscopy and

TS ranking, and the absolute best according to DFT-MBD ranking. Thus, our crystal structure prediction protocol was successful in this case. However, the statistical analysis of iodine−iodine interactions in the other 132 predicted crystal structures did not reproduce their distribution observed in other molecular crystals. In both Tinker-set and FHIaims-set, the distribution of I···I contacts over the |θ1 − θ2| parameter (allowing to distinguish the Type I, Quasi-type II and Type II interactions), is very different from the one found in CSD-set. The percentage of Type II C(sp2)−I···I−C(sp2) halogen bonds in CSD amounts to 47%, whereas it is 17% in Tinker-set and 9% in FHIaims-set. One may conclude that further development of force field potentials for the description of iodine halogen bonds in molecular crystals is necessary. The crystal packing of the best predicted structure MF33 demonstrates the highest value of χ(2)tot among five predicted noncentrosymmetric structures. We also found inverse linear correlation between NLO properties of the material and I···I halogen bond length. Thus, halogen bonds contribute to NLO properties in two different ways: by favoring noncentrosymmetric molecular packing and by polarizing the halogen atom.



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.cgd.8b00529.



Complete list of predicted virtual polymorphs with their final symmetries and lattice energies, angular distribution of type I and type II of I···I interactions, and nonlinear optical properties depending on geometric parameters (PDF)

AUTHOR INFORMATION

Corresponding Author

*Phone: 1-407-374-3783. E-mail: [email protected]. ORCID

Artëm E. Masunov: 0000-0003-4924-3380 Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This work was supported by the Ministry of Education and Science of the Russian Federation (award No. 4.1157.2017/ 4.6) and by Act 211, Government of the Russian Federation (contract no. 02.A03.21.0011) in the part of NLO properties predictions. A.E.M. acknowledges support by the “Improving of the Competitiveness” program of the National Research Nuclear University MEPhI and by the Russian Science Foundation (contract no. 14-43-00052) in the part of structure predictions. This work used the Extreme Science and Engineering Discovery Environment (XSEDE) resource Stampede2 at UT Austin through allocation TGDMR180004 and the University of Central Florida Advanced Research Computing Center (https://arcc.ist.ucf.edu).



REFERENCES

(1) Beverina, L.; Pagani, G. A. pi-Conjugated zwitterions as paradigm of donor-acceptor building blocks in organic-based materials. Acc. Chem. Res. 2014, 47 (2), 319−329. H

DOI: 10.1021/acs.cgd.8b00529 Cryst. Growth Des. XXXX, XXX, XXX−XXX

Crystal Growth & Design

Article

theoretic study of hyperpolarizability effect in diiodobutenyl-bisthioquinolinium triiodide at low temperature. J. Raman Spectrosc. 2017, 48 (11), 1411−1413. (21) Desiraju, G. R. Supramolecular synthons in crystal engineering a new organic-synthesis. Angew. Chem., Int. Ed. Engl. 1995, 34 (21), 2311−2327. (22) Nangia, A.; Desiraju, G. R. Supramolecular synthons and pattern recognition. Top. Curr. Chem. 1998, 198, 57−95. (23) Hofmann, D. W. M.; Kuleshova, L. N.; Antipin, M. Y. Supramolecular synthons and crystal structure prediction of organic compounds. Cryst. Growth Des. 2004, 4 (6), 1395−1402. (24) Yang, Y.; Xue, M.; Xiang, J. F.; Chen, C. F. Noncovalent synthesis of shape-persistent cyclic hexamers from ditopic hydrazidebased supramolecular synthons and asymmetric induction of supramolecular chirality. J. Am. Chem. Soc. 2009, 131 (35), 12657− 12663. (25) Masunov, A. E.; Tannu, A.; Dyakov, A. A.; Matveeva, A. D.; Freidzon, A. Y.; Odinokov, A. V.; Bagaturyants, A. A. First principles crystal engineering of nonlinear optical materials. I. Prototypical case of urea. J. Chem. Phys. 2017, 146 (24), 244104. (26) Reilly, A. M.; Cooper, R. I.; Adjiman, C. S.; Bhattacharya, S.; Boese, A. D.; Brandenburg, J. G.; Bygrave, P. J.; Bylsma, R.; Campbell, J. E.; Car, R.; Case, D. H.; Chadha, R.; Cole, J. C.; Cosburn, K.; Cuppen, H. M.; Curtis, F.; Day, G. M.; DiStasio, R. A.; Dzyabchenko, A.; van Eijck, B. P.; Elking, D. M.; van den Ende, J. A.; Facelli, J. C.; Ferraro, M. B.; Fusti-Molnar, L.; Gatsiou, C. A.; Gee, T. S.; de Gelder, R.; Ghiringhelli, L. M.; Goto, H.; Grimme, S.; Guo, R.; Hofmann, D. W. M.; Hoja, J.; Hylton, R. K.; Iuzzolino, L.; Jankiewicz, W.; de Jong, D. T.; Kendrick, J.; de Klerk, N. J. J.; Ko, H. Y.; Kuleshova, L. N.; Li, X. Y.; Lohani, S.; Leusen, F. J. J.; Lund, A. M.; Lv, J.; Ma, Y. M.; Marom, N.; Masunov, A. E.; McCabe, P.; McMahon, D. P.; Meekes, H.; Metz, M. P.; Misquitta, A. J.; Mohamed, S.; Monserrat, B.; Needs, R. J.; Neumann, M. A.; Nyman, J.; Obata, S.; Oberhofer, H.; Oganov, A. R.; Orendt, A. M.; Pagola, G. I.; Pantelides, C. C.; Pickard, C. J.; Podeszwa, R.; Price, L. S.; Price, S. L.; Pulido, A.; Read, M. G.; Reuter, K.; Schneider, E.; Schober, C.; Shields, G. P.; Singh, P.; Sugden, I. J.; Szalewicz, K.; Taylor, C. R.; Tkatchenko, A.; Tuckerman, M. E.; Vacarro, F.; Vasileiadis, M.; Vazquez-Mayagoitia, A.; Vogt, L.; Wang, Y. C.; Watson, R. E.; de Wijs, G. A.; Yang, J.; Zhu, Q.; Groom, C. R. Report on the sixth blind test of organic crystal structure prediction methods. Acta Crystallogr., Sect. B: Struct. Sci., Cryst. Eng. Mater. 2016, 72, 439−459. (27) Bardwell, D. A.; Adjiman, C. S.; Arnautova, Y. A.; Bartashevich, E.; Boerrigter, S. X. M.; Braun, D. E.; Cruz-Cabeza, A. J.; Day, G. M.; Della Valle, R. G.; Desiraju, G. R.; van Eijck, B. P.; Facelli, J. C.; Ferraro, M. B.; Grillo, D.; Habgood, M.; Hofmann, D. W. M.; Hofmann, F.; Jose, K. V. J.; Karamertzanis, P. G.; Kazantsev, A. V.; Kendrick, J.; Kuleshova, L. N.; Leusen, F. J. J.; Maleev, A. V.; Misquitta, A. J.; Mohamed, S.; Needs, R. J.; Neumann, M. A.; Nikylov, D.; Orendt, A. M.; Pal, R.; Pantelides, C. C.; Pickard, C. J.; Price, L. S.; Price, S. L.; Scheraga, H. A.; van de Streek, J.; Thakur, T. S.; Tiwari, S.; Venuti, E.; Zhitkov, I. K. Towards crystal structure prediction of complex organic compounds - a report on the fifth blind test. Acta Crystallogr., Sect. B: Struct. Sci. 2011, 67, 535−551. (28) Glass, C. W.; Oganov, A. R.; Hansen, N. USPEX - Evolutionary crystal structure prediction. Comput. Phys. Commun. 2006, 175 (11− 12), 713−720. (29) Dong, H. F.; Oganov, A. R.; Zhu, Q.; Qian, G. R. The phase diagram and hardness of carbon nitrides. Sci. Rep. 2015, 5, 9870. (30) Dong, H. F.; Oganov, A. R.; Wang, Q. G.; Wang, S. N.; Wang, Z. H.; Zhang, J.; Esfahani, M. M. D.; Zhou, X. F.; Wu, F. G.; Zhu, Q. Prediction of a new ground state of superhard compound B6O at ambient conditions. Sci. Rep. 2016, 6, DOI: 10.1038/srep31288. (31) Gavryushkin, P. N.; Behtenova, A.; Popov, Z. I.; Bakakin, V. V.; Likhacheva, A. Y.; Litasov, K. D.; Gavryushkin, A. Toward analysis of structural changes common for alkaline carbonates and binary compounds: prediction of high-pressure structures of Li2CO3, Na2CO3, and K2CO3. Cryst. Growth Des. 2016, 16 (10), 5612−5617.

(32) Scheiner, S. Detailed comparison of the pnicogen bond with chalcogen, halogen, and hydrogen bonds. Int. J. Quantum Chem. 2013, 113 (11), 1609−1620. (33) Scheiner, S. Sensitivity of noncovalent bonds to intermolecular separation: hydrogen, halogen, chalcogen, and pnicogen bonds. CrystEngComm 2013, 15 (16), 3119−3124. (34) Masunov, A. E.; Zorkii, P. M. Geometric characteristics of halogen-halogen intermolecular contacts in organic-crystals. Zh. Fiz. Khim. 1992, 66 (1), 60−69. (35) Bond, A. D.; Griffiths, J.; Rawson, J. M.; Hulliger, J. Inducing structural polarity using fluorinated organics: X-ray crystal structures of p-XC6F4CN (X = Cl, Br, I). Chem. Commun. 2001, 23, 2488− 2489. (36) Aakeroy, C. B.; Spartz, C. L.; Dembowski, S.; Dwyre, S.; Desper, J. A systematic structural study of halogen bonding versus hydrogen bonding within competitive supramolecular systems. IUCrJ 2015, 2, 498−510. (37) Mukherjee, A.; Tothadi, S.; Desiraju, G. R. Halogen bonds in crystal engineering: like hydrogen bonds yet different. Acc. Chem. Res. 2014, 47 (8), 2514−2524. (38) Masunov, A. E.; Grishchenko, S. I.; Zorkii, P. M. Effect of specific intermolecular interactions on crystalline-structure - biphenyl para-substituted derivatives. Zh. Fiz. Khim. 1992, 66 (1), 46−59. (39) Varughese, S. Non-covalent routes to tune the optical properties of molecular materials. J. Mater. Chem. C 2014, 2 (18), 3499−3516. (40) Marini, A.; Macchi, S.; Jurinovich, S.; Catalano, D.; Mennucci, B. Integrated NMR and computational study of push-pull NLO probes: interplay of solvent and structural effects. J. Phys. Chem. A 2011, 115 (35), 10035−10044. (41) Saccone, M.; Dichiarante, V.; Forni, A.; Goulet-Hanssens, A.; Cavallo, G.; Vapaavuori, J.; Terraneo, G.; Barrett, C. J.; Resnati, G.; Metrangolo, P.; Priimagi, A. Supramolecular hierarchy among halogen and hydrogen bond donors in light-induced surface patterning. J. Mater. Chem. C 2015, 3 (4), 759−768. (42) Chen, T. L.; Sun, Z. H.; Liu, X. T.; Wang, J. Y.; Zhou, Y. L.; Ji, C. M.; Zhang, S. Q.; Li, L. N.; Chen, Z. N.; Luo, J. H. Strong enhancement of second harmonic generation in nonlinear optical crystals: 2-amino-3-nitropyridinium halides (Cl, Br, I). J. Mater. Chem. C 2014, 2 (41), 8723−8728. (43) Virkki, M.; Tuominen, O.; Forni, A.; Saccone, M.; Metrangolo, P.; Resnati, G.; Kauranen, M.; Priimagi, A. Halogen bonding enhances nonlinear optical response in poled supramolecular polymers. J. Mater. Chem. C 2015, 3 (13), 3003−3006. (44) Desiraju, G. R.; Parthasarathy, R. The nature of halogen··· halogen interactions - are short halogen contacts due to specific attractive forces or due to close packing of nonspherical atoms. J. Am. Chem. Soc. 1989, 111 (23), 8725−8726. (45) Tothadi, S.; Joseph, S.; Desiraju, G. R. Synthon modularity in cocrystals of 4-bromobenzamide with n-alkanedicarboxylic acids: Type I and Type ll Halogen···Halogen interactions. Cryst. Growth Des. 2013, 13 (7), 3242−3254. (46) Metrangolo, P.; Resnati, G. Type II halogen···halogen contacts are halogen bonds. IUCrJ 2014, 1, 5−7. (47) Masunov, A. E.; Zorkii, P. M. Donor-acceptor nature of specific nonbonded interactions of sulfur and halogen atoms - influence on the geometry and packing of molecules. J. Struct. Chem. 1992, 33 (3), 423−435. (48) Clark, T.; Hennemann, M.; Murray, J. S.; Politzer, P. Halogen bonding: the sigma-hole. J. Mol. Model. 2007, 13 (2), 291−296. (49) Clark, T. σ-Holes. WIREs Comput. Mol. Sci. 2013, 3 (1), 13− 20. (50) Clark, T. Polarization, donor-acceptor interactions, and covalent contributions in weak interactions: a clarification. J. Mol. Model. 2017, 23 (10), 297. (51) Cavallo, G.; Metrangolo, P.; Milani, R.; Pilati, T.; Priimagi, A.; Resnati, G.; Terraneo, G. The Halogen Bond. Chem. Rev. 2016, 116 (4), 2478−2601. I

DOI: 10.1021/acs.cgd.8b00529 Cryst. Growth Des. XXXX, XXX, XXX−XXX

Crystal Growth & Design

Article

(52) Desiraju, G. R.; Ho, P. S.; Kloo, L.; Legon, A. C.; Marquardt, R.; Metrangolo, P.; Politzer, P.; Resnati, G.; Rissanen, K. Definition of the halogen bond (IUPAC Recommendations 2013). Pure Appl. Chem. 2013, 85 (8), 1711−1713. (53) Silvi, B.; Savin, A. Classification of chemical-bonds based on topological analysis of electron localization functions. Nature 1994, 371 (6499), 683−686. (54) Hunter, G. Conditional probability amplitudes in wave mechanics. Int. J. Quantum Chem. 1975, 9 (2), 237−242. (55) Vyboishchikov, S. F.; Masunov, A. E.; Streltsov, V. A.; Zorkii, P. M.; Tsirelson, V. G. Topological analysis of electron-density in chlorine crystals. Zh. Fiz. Khim. 1994, 68 (11), 2024−2028. (56) Bader, R. F. W.; Beddall, P. Virial field relationship for molecular charge distributions and the spatial partitioning of molecular properties. J. Chem. Phys. 1972, 56, 3320−3329. (57) Bartashevich, E. V.; Tsirelson, V. G. Atomic dipole polarization in charge-transfer complexes with halogen bonding. Phys. Chem. Chem. Phys. 2013, 15 (7), 2530−2538. (58) Bartashevich, E.; Troitskaya, E.; Pendas, A. M.; Tsirelson, V. Understanding the bifurcated halogen bonding N···Hal···N in bidentate diazaheterocyclic compounds. Comput. Theor. Chem. 2015, 1053, 229−237. (59) Day, G. M.; Price, S. L. A nonempirical anisotropic atom-atom model potential for chlorobenzene crystals. J. Am. Chem. Soc. 2003, 125 (52), 16434−16443. (60) Titov, O. I.; Shulga, D. A.; Palyulin, V. A.; Zefirov, N. S. Description of halogen bonding on the basis of multicenter multipole expansion. Dokl. Chem. 2013, 450, 139−143. (61) Mu, X.; Wang, Q.; Wang, L.-P.; Fried, S. D.; Piquemal, J.-P.; Dalby, K. N.; Ren, P. Modeling organochlorine compounds and the σhole effect using a polarizable multipole force field. J. Phys. Chem. B 2014, 118 (24), 6456−6465. (62) Saha, B. K.; Nangia, A.; Nicoud, J.-F. Using Halogen···Halogen interactions to direct noncentrosymmetric crystal packing in dipolar organic molecules. Cryst. Growth Des. 2006, 6 (6), 1278−1281. (63) Bartashevich, E. V.; Yushina, I. D.; Stash, A. I.; Tsirelson, V. G. Halogen bonding and other iodine interactions in crystals of dihydrothiazolo(oxazino)quinolinium oligoiodides from the electron-density viewpoint. Cryst. Growth Des. 2014, 14 (11), 5674−5684. (64) Bartashevich, E. V.; Grigoreva, E. A.; Yushina, I. D.; Bulatova, L. M.; Tsirelson, V. G. Modern level for properties prediction of iodine-containing organic compounds: the halogen bonds formed by iodine. Russ. Chem. Bull. 2017, 66 (8), 1345−1356. (65) Wu, J. C.; Chattree, G.; Ren, P. Y. Automation of AMOEBA polarizable force field parameterization for small molecules. Theor. Chem. Acc. 2012, 131 (3), 1138. (66) Towns, J.; Cockerill, T.; Dahan, M.; Foster, I.; Gaither, K.; Grimshaw, A.; Hazlewood, V.; Lathrop, S.; Lifka, D.; Peterson, G. D.; Roskies, R.; Scott, J. R.; Wilkens-Diehr, N. XSEDE: accelerating scientific discovery. Comput. Sci. Eng. 2014, 16 (5), 62−74. (67) 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.; Li, X.; Caricato, M.; Marenich, A.; Bloino, J.; Janesko, B. G.; Gomperts, R.; Mennucci, B.; Hratchian, H. P.; Ortiz, J. V.; Izmaylov, A. F.; Sonnenberg, J. L.; Williams-Young, D.; Ding, F.; Lipparini, F.; Egidi, F.; Goings, J.; Peng, B.; Petrone, A.; Henderson, T.; Ranasinghe, D.; Zakrzewski, V. G.; Gao, J.; Rega, N.; Zheng, G.; Liang, W.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Throssell, K.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Keith, T.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Millam, J. M.; Klene, M.; Adamo, C.; Cammi, R.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Farkas, O.; Foresman, J. B.; Fox, D. J. Gaussian 09, Rev. D.01; Gaussian, Inc.: Wallingford, CT, 2009.

(68) Ren, P. Y.; Ponder, J. W. Polarizable atomic multipole water model for molecular mechanics simulation. J. Phys. Chem. B 2003, 107 (24), 5933−5947. (69) Ren, P. Y.; Wu, C. J.; Ponder, J. W. Polarizable atomic multipole-based molecular mechanics for organic molecules. J. Chem. Theory Comput. 2011, 7 (10), 3143−3161. (70) Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. Selfconsistent molecular-orbital methods. 20. Basis set for correlated wave-functions. J. Chem. Phys. 1980, 72 (1), 650−654. (71) Glukhovtsev, M. N.; Pross, A.; McGrath, M. P.; Radom, L. Extension of Gaussian-2 (G2) theory to bromine-containing and iodine-containing molecules - use of effective core potentials. J. Chem. Phys. 1995, 103 (5), 1878−1885. (72) Stone, A. J. Distributed multipole analysis: Stability for large basis sets. J. Chem. Theory Comput. 2005, 1 (6), 1128−1132. (73) Zhu, Q.; Oganov, A. R.; Glass, C. W.; Stokes, H. T. Constrained evolutionary algorithm for structure prediction of molecular crystals: methodology and applications. Acta Crystallogr., Sect. B: Struct. Sci. 2012, 68, 215−226. (74) Belsky, V. K.; Zorkaya, O. N.; Zorky, P. M. Structural classes and space-groups of organic homomolecular crystals - new statisticaldata. Acta Crystallogr., Sect. A: Found. Crystallogr. 1995, 51, 473−481. (75) Stokes, H. T.; Hatch, D. M. FINDSYM: Program for identifying the space group symmetry of a crystal. J. Appl. Crystallogr. 2005, 38, 237−238. (76) Macrae, C. F.; Bruno, I. J.; Chisholm, J. A.; Edgington, P. R.; McCabe, P.; Pidcock, E.; Rodriguez-Monge, L.; Taylor, R.; van de Streek, J.; Wood, P. A. Mercury CSD 2.0-new features for the visualization and investigation of crystal structures. J. Appl. Crystallogr. 2008, 41, 466−470. (77) Blum, V.; Gehrke, R.; Hanke, F.; Havu, P.; Havu, V.; Ren, X. G.; Reuter, K.; Scheffler, M. Ab initio molecular simulations with numeric atom-centered orbitals. Comput. Phys. Commun. 2009, 180 (11), 2175−2196. (78) Tkatchenko, A.; Scheffler, M. Accurate molecular van der waals interactions from ground-state electron density and free-atom reference data. Phys. Rev. Lett. 2009, 102 (7), 073005. (79) Tkatchenko, A.; DiStasio, R. A.; Car, R.; Scheffler, M. Accurate and efficient method for many-body van der waals interactions. Phys. Rev. Lett. 2012, 108 (23), 236402. (80) Marom, N.; DiStasio, R. A.; Atalla, V.; Levchenko, S.; Reilly, A. M.; Chelikowsky, J. R.; Leiserowitz, L.; Tkatchenko, A. Many-body dispersion interactions in molecular crystal polymorphism. Angew. Chem., Int. Ed. 2013, 52 (26), 6629−6632. (81) Groom, C. R.; Bruno, I. J.; Lightfoot, M. P.; Ward, S. C. The Cambridge Structural Database. Acta Crystallogr., Sect. B: Struct. Sci., Cryst. Eng. Mater. 2016, 72, 171−179. (82) Dovesi, R.; Orlando, R.; Erba, A.; Zicovich-Wilson, C. M.; Civalleri, B.; Casassa, S.; Maschio, L.; Ferrabone, M.; De La Pierre, M.; D’Arco, P.; Noel, Y.; Causa, M.; Rerat, M.; Kirtman, B. CRYSTAL14: A Program for the ab initio investigation of crystalline solids. Int. J. Quantum Chem. 2014, 114 (19), 1287−1317. (83) Grimme, S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction. J. Comput. Chem. 2006, 27 (15), 1787−1799. (84) Barros, C. L.; de Oliveira, P. J. P.; Jorge, F. E.; Canal Neto, A.; Campos, M. Gaussian basis set of double zeta quality for atoms Rb through Xe: application in non-relativistic and relativistic calculations of atomic and molecular properties. Mol. Phys. 2010, 108 (15), 1965− 1972. (85) Gatti, C.; Saunders, V. R.; Roetti, C. Crystal-field effects on the topological properties of the electron-density in molecular-crystals the case of urea. J. Chem. Phys. 1994, 101 (12), 10686−10696. (86) Maschio, L.; Kirtman, B.; Rerat, M.; Orlando, R.; Dovesi, R. Ab initio analytical Raman intensities for periodic systems through a coupled perturbed Hartree-Fock/Kohn-Sham method in an atomic orbital basis. I. Theory. J. Chem. Phys. 2013, 139 (16), 164101. (87) Keshari, V.; Wijekoon, W. M. K. P.; Prasad, P. N.; Karna, S. P. Hyperpolarizabilities of organic-molecules - ab-initio time-dependent J

DOI: 10.1021/acs.cgd.8b00529 Cryst. Growth Des. XXXX, XXX, XXX−XXX

Crystal Growth & Design

Article

coupled perturbed hartree-fock-roothaan studies of basic heterocyclic structures. J. Phys. Chem. 1995, 99 (22), 9045−9050. (88) Kurtz, S. K.; Perry, T. T. A powder technique for the evaluation of nonlinear optical materials. J. Appl. Phys. 1968, 39, 3798. (89) Pereira Silva, P. S.; Pereira Gonçalves, M. A.; Ramos Silva, M.; Paixão, J. A. Structural and nonlinear optical studies of a salt with an octupolar chromophore: Guanidinium cyclopropanecarboxylate. Spectrochim. Acta, Part A 2017, 172, 156−162. (90) Gatti, C.; Casassa, S. TOPOND14 User Manual. http://www. crystal.unito.it/topond/topond.pdf (2016). (91) Bartashevich, E.; Yushina, I.; Kropotina, K.; Muhitdinova, S.; Tsirelson, V. Testing the tools for revealing and characterizing the iodine-iodine halogen bond in crystals. Acta Crystallogr., Sect. B: Struct. Sci., Cryst. Eng. Mater. 2017, 73, 217−226. (92) Bol’shakov, O. I.; Yushina, I. D.; Bartashevich, E. V.; Nelyubina, Y. V.; Aysin, R. R.; Rakitin, O. A. Asymmetric triiodide-diiodine interactions in the crystal of (Z)-4-chloro-5-((2-((4-chloro-5H-1,2,3dithiazol-5-ylidene)amino)phenyl)amino)-1,2,3-dithiazol-1-ium oligoiodide. Struct. Chem. 2017, 28 (6), 1927−1934. (93) Bader, R. F. W. Atoms in Molecules: A Quantum Theory. Clarendon Press: Oxford, 1990.

K

DOI: 10.1021/acs.cgd.8b00529 Cryst. Growth Des. XXXX, XXX, XXX−XXX