Intermolecular Interactions in Ionic Crystals of Nucleobase Chlorides—Combining Topological Analysis of Electron Densities with Energies of Electrostatic Interactions

Abstract

Understanding intermolecular interactions in crystals of molecular ions continues to be difficult. On the one hand, the analysis of interactions from the point of view of formal charges of molecules, similarly as it is commonly done for inorganic ionic crystals, should be performed. On the other hand, when various functional groups are present in the crystal, it becomes natural to look at the interactions from the point of view of hydrogen bonding, π…π stacking and many other kinds of non-covalent atom–atom bonding. Often, these two approaches seem to lead to conflicting conclusions. On the basis of experimental charge densities of cytosinium chloride, adeninium chloride hemihydrate, and guanine dichloride crystals, with the help of theoretical simulations, we have deeply analysed intermolecular interactions among protonated nucleobases, chloride anions and water molecules. Here, in the second paper of the series of the two (Kumar et al., 2018, IUCrJ 5, 449–469), we focus on applying the above two approaches to the large set of dimers identified in analysed crystals. To understand electrostatic interactions, we analysed electrostatic interaction energies (E_es) computed directly from molecular charge densities and contrasted them with energies computed only from net molecular charges, or from a sum of electric multipolar moments, to find the charge penetration contribution to E_es. To characterize non-covalent interactions we performed topological analyses of crystal electron densities and estimated their interaction energies (E_EML) from properties of intermolecular bond critical points. We show that the overall crystal architecture of the studied compounds is governed by the tight packing principle and strong electrostatic attractions and repulsions between ions. Many ions are oriented to each other in a way to strengthen attractive electrostatic interactions or weaken strong repulsion, but not all of them. Numerous bond critical points and bond paths were found between ions, including nucleobase cations despite their overall repulsive interactions. It is clear there is no correlation between E_EML and E_es. However, strong relation between E_EML and the charge penetration component of E_es is observed. The relation holds regardless of interaction types or whether or not interacting molecules bear the same or opposite charges. Thus, a charge density-based approach for computing intermolecular interaction energies and the atom–atom approach to analyse non-covalent interactions do complement each other, even in ionic systems.

1. Introduction

Intermolecular interactions play a crucial role in many areas of chemistry and biology, e.g., synthetic chemistry, crystal engineering, protein folding, gene expression and drug discovery. They are responsible for molecular recognition and supramolecular assembly of macromolecular (e.g., secondary and higher structures of proteins or nucleic acids) or small molecular (e.g., crystals) systems. A detailed understanding of these interactions has outstanding importance in rationalizing structure–function and structure–stability relationships observed in these fields.

In the field of crystal engineering, understanding intermolecular interactions is essential in the context of crystal packing in order to use it to design new solids with desired physical and chemical properties [1]. In 1995 Professor Desiraju introduced the term supramolecular synthon. Supramolecular synthons are essentially structural units within supramolecules which can be formed (assembled) by known or conceivable operations involving non-covalent interactions [2]. These are kinetically defined structural units (building blocks) which arrange themselves in favour of strong and directional interactions (atom–atom interaction) to form a crystal. Later, following the idea of synthons, Dunitz and Gavezzotti performed PIXEL energy calculations on various systems and confirmed that synthons correspond to stable crystal building blocks [3]. However, the authors also infer that this approach is not straightforward since crystal structures can be broken down in many ways. In addition many, almost equi-energetic crystal structures (e.g., in polymorphs) are possible for a given compound, which further complicates the situation. In particular, for crystals involving weaker interactions, many calculated energy surfaces are quite flat, so the exact geometry of the synthon may be variable and one-coordinate descriptions in terms of a single atom–atom bond is inadequate. Hence, one must look into interactions between whole molecular charge densities, rather than pairwise atom–atom interactions only, especially at short intermolecular distances or in the presence of highly charged groups [3,4,5,6,7,8,9].

The situation becomes even more complicated in ionic crystals where anion–anion (A⁻...A⁻), cation–cation (C⁺…C⁺) and anion–cation (A⁻…C⁺) interactions are the major types of non-covalent interactions involved in crystal stabilization. According to classical molecular mechanics, ionic crystals are characterized by strong, stabilizing or destabilizing interactions depending on the attractive or repulsive nature of the forces acting between the interacting molecules. Opposite charge pairing (A⁻…C⁺) is easily understood and one of the most widespread ideas in chemistry, ion pairs or ionic interactions are often responsible for chemical reactions or, more particularly, are important in designing new materials. Over the past few years, considerable effort has been devoted to comprehending the existence of cation–cation (C⁺…C⁺) and anion-anion (A⁻...A⁻) aggregation by hydrogen bond interactions in the solid state despite ionic repulsion [10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26]. To elaborate on this, Braga and co-workers [17,20,25,27] studied the ionic crystal structure of hydrogen oxalate, KHC₂O₄, alkali hydrogen squarate salts and alkali metal croconate salts, and introduced a concept called Coulombic compression ”arising from attractive next-neighbour A⁻…C+ interactions, which largely overcompensate the combined effect of next-neighbour A⁻...A⁻ and C⁺…C⁺ repulsions and allows deeper penetration into the repulsive walls of neighbouring atoms”. Hence, ions aggregate with high complexity and equilibrium is reached at the balancing point between attractive and repulsive forces, with net stabilization energy resulting from global (free) energy minimums in the solid state [9,25,27]. Various computational concepts were introduced to deeper understand attraction–repulsion balance among molecules in the solid state [28,29].

The present paper is the second in a series of two. The series is dedicated to the studies of intermolecular interactions of the protonated cytosine, adenine and guanine nucleobases with themselves, chloride anions or water molecules observed in their crystalline state by X-ray diffraction (Scheme 1) complemented with quantum mechanics calculations. Our aim is to qualitatively and quantitatively characterise these interactions in order to understand the interplay between the molecular and crystal architecture in ionic systems. In the first paper [30] of the series we provided detailed structural analysis of molecular motifs found in the studied crystals to give a basis for further charge density and energy investigations. We discovered that molecules in the studied crystals are not fully ionized and some dimers of single protonated bases, namely cytosinium–cytosinium trans sugar/sugar edge pairs and adeninium–adeninium trans Hoogsteen/Hoogsteen edge pairs, exhibited attractive electrostatic interactions or unusually low repulsion despite identical molecular charges. We concluded that these two pairs “are metastable as a result of strong hydrogen bonding between bases which overcompensate overall cation–cation repulsion, the latter being weakened due to charge transfer and molecular charge density polarization”.

Here we continue a detailed analysis of electrostatic intermolecular interactions for a larger variety of molecular dimers identified in the studied crystals and provide a comprehensive topological analysis of observed atom–atom intermolecular bonding. We show how short (at single atom, functional group and whole molecule levels) and long range (with the second nearest neighbour molecules and further away) interactions contribute to the stability of the studied ionic crystals. We aim to find a consistent view of interactions in ionic molecular crystals from both approaches: topology analysis of electron densities focused on atom–atom contacts and interactions energies computed for entire interacting molecules. All of these to deeper understand how molecular ions pack together to form larger supramolecular assemblies, like crystal or biomacromolecular complexes.

2. Materials and Methods

2.1. Crystal Charge Density Models

The following charge density models of cytosinum chloride (CC), adeninium chloride hemihydrate (ACH) and guaninium dichloride (GDC) crystals were taken from [30] and continued to be analysed in this paper:

(a) Experimental charge densities (Exp.) in the Stewart–Hansen–Coppens [31,32] multipole model (MM) representation with the Cl⁻¹ ion scattering radial function and ion configuration;

(b) Charge densities reconstructed with University at Buffalo Databank (UBDB) in the Stewart–Hansen–Coppens multipole model representation parameterised on experimental geometries using the LSDB program [33] and the UBDB2011 version of the databank [34];

(c) Theoretical electron densities computed directly from periodic wave functions obtained from periodic quantum mechanical geometry optimization (Theo. opt.) or from single-point calculations at experimental geometry (Theo. Exp.) done with the Crystal14 package [35,36] at the DFT-B3LYP/pVDZ level of theory [37,38,39,40];

(d) Theoretical charge densities in Stewart–Hansen–Coppens multipole model representation (Theo. opt. MM; Theo. exp. MM) with the Cl⁻¹ ion scattering radial function and ion configuration obtained from multipolar refinement against theoretical structure factors computed from (c).

2.2. QTAIM Analysis

For models in the multipole representation (Exp., UBDB, Theo. opt. MM and Theo. exp. MM) bond critical points (BCPs) were found and characterized using the TOPXD module of the XD2016 package [41]. For exact (not approximated by the multipolar model) theoretical crystal electron densities (Theo. opt.) topological analysis was done using TOPOND14 [42]. To prepare Figure 1 and Figures S1–S3 in SI the MoPro package [43,44] was used. Topology analysis done in MoPro led to the same results as the one from XD2016.

2.3. Intermolecular Interaction Energies and Electrostatic Contributions to Them

The following interaction energies for isolated dimers with geometry as found in studied crystals were taken from [30] and continued to be analysed in this paper:

(a) The total intermolecular interaction energies $E_{tot}$ and exact electrostatic energies $E_{es}$ (originally termed $E_{Pol}^1$) computed at experimental geometry within the DFT-based symmetry-adapted perturbation theory (DFT-SAPT) [45,46,47] in MOLPRO2012.1 [48] at the PBE0AC/aug-cc-pVTZ level of theory [40,49,50,51];

(b) The exact electrostatic energies $E_{es}$ computed at experimental geometry within the semi-classical density sum methodology (the PIXEL method) [52] from MP2/6-31G** [53] molecular wave functions using the PIXELc module [54,55,56,57] of the Coulomb–London–Pauli (CLP) program [52]. The program allows to compute in automatic mode energies between the central molecule and any other molecule surrounding it up to large distance, what is important for ionic systems and is not easy to be done in other programs like CrystalExplorer17 [58], for example.

(c) The electrostatic energies of various kinds, as explained below, computed from multipolar models of experimental (Exp.), UBDB and periodic theoretical charge densities (Theo. opt. MM) at respective geometries, using the XDPROP module of the XD2006 package [41]. Electrostatic interaction contributions were computed applying several levels of approximation. The exact electrostatic energies ($E_{es}^{ }$) and the electrostatic energies from molecular monopole moments (point charges,$E_{Coul}^{}$) were taken directly from [30]. In addition, the electrostatic energies from multipole expansion were computed from: atomic multipole moments ($E_{mtp}^{atom}$) using the aMM option and from molecular multipole moments ($E_{mtp}^{mol}$) using the mMM option.

3. Results and Discussion

Firstly, detailed analysis of intermolecular interactions in studied crystals was performed from the geometrical point of view on the basis of interatomic distances, angles and van der Waals radii. The results are given in the first article from the series; see [30].

3.1. Topological Analysis of Electron Density

Topological analysis of experimental electron density, reported here, enriched our knowledge of intermolecular interactions in the crystal structures. In general, considerably more bond paths were identified using topological analysis than one could expect from simple geometrical analysis based on van der Waals radii; see Figure 1,

Table 1. Selected QTAIM parameters of intermolecular interactions at bond critical points (BCPs) detected in the experimental electron densities of the CC, ACH and GDC crystals. R_ab—distance (Å) between interacting atoms; d1_BPL and d2_BPL—bond path lengths (Å) from the first or the second atom to the BCP, respectively; ρ(r)_BCP—electron density (eÅ⁻³) at BCP; ∇²ρ(r)_BCP—Laplacian of electron density (eÅ⁻⁵) at BCP; ε—ellipticity at BCP, E_EML—interaction energy (kcal mol⁻¹) estimated from Espinosa-Molins-Lecomte approach. For the symmetry operations required to build particular dimers see Table S2 in SI. In the Interaction type column, HB denotes a hydrogen bond. In bold atom–atom interactions and dimers identified on the basis of the van der Waals radii, see Table 3 in [30].

Dimer	No. BCP	Interacting Atoms	Interaction Type	Rab	d1BPL	d2BPL	ρ(r)BCP	∇2ρ(r)BCP	ε	EEML
Cytosinium chloride (CC)
AA1	CP1	H1...O2	HB	1.776	0.640	1.136	0.307	1.64	0.07	−12.2
AA2	CP2	C4...C2	π...π	3.312	1.668	1.752	0.047	0.43	0.72	−0.9
AA2	CP3	N3...N3	π...π	3.312	1.656	1.656	0.039	0.52	14.39	−0.9
AA3	CP4	H6...O2	HB	2.327	0.955	1.392	0.082	1.11	0.36	−2.4
AA4	CP5	N7...C5	π...π	3.433	1.688	1.797	0.029	0.39	0.65	−0.6
AA5	CP6	N7...N1	π...π	3.765	1.893	1.905	0.011	0.17	1.54	−0.2
AA6	CP7	H6...C5	CH... π	3.007	1.285	1.802	0.030	0.28	0.20	−0.5
AB1	CP8	H7B...Cl1	HB	2.362	0.867	1.508	0.138	0.86	0.05	−3.7
AB1	CP9	H5...Cl1	HB	2.545	0.924	1.633	0.084	1.01	0.29	−2.3
AB2	CP10	N7...Cl1	NH2...Cl-	3.308	1.575	1.738	0.051	0.74	0.63	−1.3
AB3	CP11	C2...Cl1	π...Cl−	3.340	1.578	1.767	0.053	0.58	0.32	−1.2
AB4	CP12	H3...Cl1	HB	2.022	0.690	1.337	0.304	1.46	0.05	−11.9
AB4	CP13	H7A...Cl1	HB	2.619	0.973	1.699	0.063	0.80	0.58	−1.6
AB5	CP14	C6...Cl1	π...Cl−	3.729	1.932	2.003	0.024	0.36	0.68	−0.5
AB6	CP15	C6...Cl1	π...Cl−	3.691	1.840	1.946	0.024	0.29	0.40	−0.5
Adeninium chloride hemihydrate (ACH)
AA1	CP1	H10A...N7	HB	1.937	0.708	1.235	0.225	1.67	0.03	−8.0
AA2	CP2	H2...N3	HB	2.734	1.358	1.562	0.048	0.72	0.20	−1.2
AA2	CP3	N3...N3		3.150	1.585	1.585	0.040	0.67	0.53	−1.1
AA3	−
AA4	CP4	N3...N7	π...π	3.217	1.609	1.619	0.043	0.56	0.81	−1.0
AA4	CP5	C2...N10	π...π	3.438	1.767	1.917	0.030	0.37	3.05	−0.6
AA5	CP6	N7…N10	π...π	3.340	1.631	1.980	0.045	0.52	1.15	−1.0
AB1	CP7	H10B...Cl1	HB	2.255	0.785	1.480	0.144	1.26	0.06	−4.3
AB2	CP8	H8...Cl1	HB	2.606	0.936	1.683	0.066	0.65	0.15	−1.5
AB3	CP9	H9...Cl1	HB	2.109	0.738	1.378	0.221	1.04	0.00	−7.2
AB4	CP10	H2...Cl1	HB	2.699	0.999	1.706	0.054	0.65	0.04	−1.3
AB5	CP11	N1...Cl1	π...Cl−	3.569	1.724	1.866	0.027	0.38	0.67	−0.6
AB6	CP12	C8...Cl1	π...Cl−	3.950	1.952	2.083	0.012	0.15	0.42	−0.2
AB7	CP13	N9...Cl1	π...Cl−	3.826	1.812	2.035	0.014	0.22	3.09	−0,3
AW1	CP14	H1...O1	HB	1.834	0.678	1.159	0.268	2.01	0.02	−10.5
BW1	CP15	H1A...Cl1	HB	2.118	0.704	1.418	0.217	0.96	0.02	−6.9
BW2	CP16	O1...Cl1		3.608	1.680	1.932	0.021	0.32	0.67	−0.5
BB1	CP17	Cl1...Cl1	Cl−...Cl−	4.262	2.132	2.132	0.010	0.14	0.09	−0.2
Guaninium dichloride (GDC)
AA1	CP1	H8...O10	HB	2.246	1.179	1.304	0.083	1.69	0.08	−3.1
AA1	CP2	N9...O10		2.759	1.725	1.324	0.080	1.52	0.65	−2.8
AA2	CP3	O10...N11	π...π	3.418	1.731	1.737	0.023	0.33	0.72	−0.5
AB1	CP4	H11B...Cl1	HB	2.084	0.676	1.413	0.193	1.54	0.07	−6.5
AB1	CP5	H3...Cl1	HB	2.274	0.821	1.467	0.151	0.98	0.03	−4.3
AB2	CP6	H1...Cl1	HB	2.036	0.664	1.373	0.231	2.13	0.02	−8.8
AB3	CP7	H8...Cl1	HB	2.826	1.083	1.758	0.047	0.27	0.05	−0.6
AB4	CP8	H3...Cl2	HB	2.996	1.246	1.851	0.030	0.35	1.27	−0.6
AB4	CP9	H9...Cl2	HB	2.204	0.765	1.445	0.162	1.37	0.02	−5.1
AB5	CP10	H7...Cl2	HB	2.076	0.712	1.367	0.262	1.27	0.01	−9.4
AB6	CP11	H11A...Cl2	HB	2.274	0.805	1.478	0.145	1.11	0.04	−4.2
AB7	CP12	C8...Cl1	π...Cl−	3.255	1.517	1.747	0.049	0.60	0.75	−1.1
AB8	CP13	C2...Cl2	π...Cl−	3.424	1.643	1.826	0.038	0.47	2.23	−0.8
AB8	CP14	C4...Cl2	π...Cl−	3.546	1.951	1.833	0.036	0.43	3.26	−0.8
AB8	CP15	N1...Cl2	π...Cl−	3.483	1.684	1.815	0.037	0.51	4.03	−0.9
BB1	CP16	Cl1...Cl2	Cl−...Cl−	3.692	1.846	1.848	0.031	0.40	0.73	−0.7

and Figures S1–S3 in SI and compare with Table 3 in [30]. Nonetheless, bond paths were not always found for atom pairs between which the distance is shorter than the sum of van der Waals radii. The abundance of bond paths indicate that a larger number of molecules directly interact with each other and more types of dimers representing local interactions should be taken into account in order to understand the architecture of studied crystal structures. The analysis of electron density topology confirms what was already assumed from the Hirshfeld surface analysis presented in [30].

Bond paths (BPs) and bond critical points (BCPs) were found (see HB in Table 1) for almost all hydrogen bonds previously identified on the basis of geometry (see HB in Table 3 in [30]), regardless if they are formed between cations and anions or between two cations. The range of values of BCP properties are similar to the ones for neutral base pairs [59,60]. One new hydrogen bond was recognized in the GDC structure, i.e., N3–H3...Cl2 in dimer AB4. The N7–H7B...Cl1 hydrogen bond in dimer AB2 of the CC structure was not confirmed by topological analysis, instead a bond path between N7 and Cl1 was found. The N9–H9...O10 hydrogen bond in dimer AA1 of the GDC structure also seems not to be present according to topological analysis. A search along the bond path originating from BCP found, in between H9 and O10, atoms identified the N9 atom as an attractor, not the H9, but the shape of the path is unusual; see Figure 1. To have more evidence regarding hydrogen bonding in the studied structures we performed topological analysis for electron densities directly obtained from periodic quantum mechanical calculations (not through reciprocal-space fitting of MM). The periodic calculations were done either with geometry optimization (Theo. opt.), or for experimental geometry (Theo. exp.). Theoretical results confirmed all the above findings with the exception of the N9–H9...O10 hydrogen bond. For that particular hydrogen bond the bond path between H9 and O10 was found as expected (Table S1 in SI). To further understand the source of the observed discrepancy between experiment and theory, we performed topological analyses of theoretical charge densities in multipole model representation (Theo. opt. MM and Theo. exp. MM). In case of Theo. opt. MM the H9…O10 bond path was found, but for Theo. exp. MM, similarly like for experimental density, the bond path was located between N9 and O10. Apparently, experimental charge density suffers somewhat from inadequate accuracy in the hydrogen atom location followed by difficulties in electron density modelling through the multipole model.

Judging from values of $\rho \left(r\right)$ at BCPs being in the range of 0.31 eÅ⁻³ to 0.03 eÅ⁻³, hydrogen bonds, including those in which carbon atoms act as proton donors, are the strongest interaction types among all intermolecular interactions found in the studied structures. The experimental values of $\rho {\left(r\right)}_{BCP}$ agree with the theoretical ones obtained for electron densities from periodic quantum mechanical calculations combined with structure relaxation (geometry optimization) (Figure S4 and Table S1 in SI). The difference does not exceed 0.04 eÅ⁻³ for the strongest interactions, and is not significantly larger than the experimental standard deviation errors being in the range of 0.01 eÅ⁻³ to 0.03 eÅ⁻³. It is interesting to note that values of $\rho {\left(r\right)}_{BCP}$ for the UBDB electron densities are systematically lower when compared to experimental results. The UBDB model does not take into account electron density response to intermolecular interactions because it is based on non-interacting model molecules. In addition, strong disagreement between the experimental and any other theoretical method (UBDB in MM representation or periodic DFT without MM) for values of ${\nabla }^2\rho {\left(r\right)}_{BCP}$ of short NH/OH…acceptor contacts is observed (Figure S4 and Table S1 in SI).

Despite anomalous values of ${\nabla }^2\rho {\left(r\right)}_{BCP}$ from the experimental density, it is tempting to compute the interaction energies ($E_{EML}$) estimated from $\rho {\left(r\right)}_{BCP}$ and ${\nabla }^2\rho {\left(r\right)}_{BCP}$ on the basis of Espinosa–Molins–Lecomte (EML) approximation [61] and Abramov expression [62]. It appears that hydrogen bonds of the O/N–H...Cl⁻ type identified as particularly strong on the basis of geometric features have indeed the largest negative values of $E_{EML}$ energies, in the range of −7 to −12 kcal mol⁻¹. Interestingly, similarly negative values of $E_{EML}$ interaction energies are found for some hydrogen bonds linking two nucleobase cations, like N1–H1...O2 in the dimer AA1 of CC or N10–H10A...N7 in the dimer AA1 of ACH, for which overall cation...cation interaction energy should be positive due to repulsion between like charges. It is to be noted that, despite noticeable difference in values of ${\nabla }^2\rho {\left(r\right)}_{BCP}$, $E_{EML}$ energies based on experimental data differ on average only by 0.8 kcal mol⁻¹ from $E_{EML}$ energies approximated from the exact potential energy densities derived from periodic quantum mechanical calculations combined with geometry optimisation (Table S1 in SI).

In addition to X–H…Cl⁻ hydrogen bonding, another type of direct interaction between nucleobase cations and chloride anions, namely π…Cl⁻, is corroborated by topological analysis. The analysis indicates that there are more such interactions than previously identified on the basis of geometry. Thus, in the CC structure, in addition to dimer AB3, two other types of dimers with π…Cl⁻ interactions are found (AB5 and AB6, Table 1). In ACH, in which no π...Cl⁻ was previously noted, three such interactions are identified (dimers AB5, AB6 and AB7). In GDC, π...Cl⁻ interaction involving the Cl1 anion is confirmed (dimer AB7), whereas in dimer AB8, bond paths from Cl2 to two more atoms from guaninium pyrimidine ring, beside C2, are located.

Next, topological analysis supports the statement that there are stacking interactions between nucleobase cations in the CC and ACH structures and shows that even in the GDC two such interactions could be identified; Table 1. In addition to previously recognized dimers AA2 in CC with cations directly interacting mainly through Watson–Crick/Watson–Crick edges, the following new dimers are identified: with cations interacting through Hoogsteen/Hoogsteen edges (AA4 in CC, AA5 in ACH), Watson–Crick/Hoogsteen edges (AA5 in CC), Hoogsteen/sugar edges (AA4 in ACH), Watson-Crick/sugar edges (AA2 in GDC) and sugar/sugar edges (AA3 in GDC). On the other hand, the presence of stacking interactions in dimers AA3 in ACH was not confirmed by topology analysis.

In addition to stacking interactions, topology analysis indicate there is one interaction of C–H...π type in the studied structures, i.e., C6–H6...C5 (dimer AA6) in CC. Another interesting type of interaction, to our knowledge not discussed in literature, found thanks to topology analysis, is a Cl⁻...O interaction in which Cl⁻ approaches the oxygen atom along the O–H direction but from the opposite side of oxygen respect to the hydrogen atom (dimer BW2 in ACH). And lastly, in ACH and GDC structures, bond paths between two chloride anions are found (dimers BB1).

Many contacts between particular chemical elements found in the Hirshfeld surface analysis [30] and not seen in the van der Waals radii based analysis are identified by the topological analysis as being sites of atom–atom interactions. But not all of them. Topological analysis does not confirm, for example, that there are important C...H interactions in the other two crystal structures, ACH and GDC. What is more crucial is the fact that, despite the abundance of H...H contacts found using Hirshfeld surface analysis, there is no a single bond path found between hydrogen atoms in any of the studied structures. To sum up, not all atom–atom contacts can be associated with atom–atom interactions according to topological analysis.

3.2. Interaction Energies for Dimers

Structural and topological analysis gave us some insights to pairwise atom–atom interactions. However to have a more global understanding of intermolecular bonding, we analysed energies of intermolecular interactions between whole molecules in dimers extracted from crystal structures,

Table 2. Total interaction energies $E_{tot}$ (kcal mol⁻¹) according to the DFT-SAPT method and electrostatic contribution (kcal mol⁻¹) to the total interaction energy for selected dimers extracted from the CC, ACH and GDC crystal structures according to different methods. $E_{\mathit{es}}$ is the electrostatic energy computed from exact integration of charge densities, $E_{mtp}$ is electrostatic energies computed from the multipole moments of the atoms ($E_{mtp}^{atom}$ ) or of the molecules ($E_{mtp}^{mol}$ ) or only from the monopole moments (overall charges) of the molecules ($E_{Coul}$ ). For the symmetry operations required to build particular dimers see Table S2 in SI. *—data taken from Table 3 in [30].

Dimer	Centre of Mass Distance (Å)	Experimental				UBDB				DFT–SAPT
		${\mathit{E}}_{\mathit{e}\mathit{s}}\ast$	${\mathit{E}}_{\mathit{m}\mathit{t}\mathit{p}}^{\mathit{a}\mathit{t}\mathit{o}\mathit{m}}$	${\mathit{E}}_{\mathit{m}\mathit{t}\mathit{p}}^{\mathit{m}\mathit{o}\mathit{l}}$	${\mathit{E}}_{\mathit{C}\mathit{o}\mathit{u}\mathit{l}}\ast$	${\mathit{E}}_{\mathit{e}\mathit{s}}\ast$	${\mathit{E}}_{\mathit{m}\mathit{t}\mathit{p}}^{\mathit{a}\mathit{t}\mathit{o}\mathit{m}}$	${\mathit{E}}_{\mathit{m}\mathit{t}\mathit{p}}^{\mathit{m}\mathit{o}\mathit{l}}$	${\mathit{E}}_{\mathit{C}\mathit{o}\mathit{u}\mathit{l}}\ast$	${\mathit{E}}_{\mathit{e}\mathit{s}}\ast$	${\mathit{E}}_{\mathit{t}\mathit{o}\mathit{t}}\ast$
Cytosinium Chloride (CC)
AA1	6.060	−0.8	17.9	18.2	42.3	14.0	24.6	30.7	54.8	16.6	20.0
AA2	3.455	42.8	45.6	46.5	74.1	60.8	63.9	108.0	96.1	61.9	55.4
AA3	6.739	27.8	28.9	29.0	38.0	40.0	41.0	40.9	49.3	39.8	38.1
AA4	5.385	57.7	58.4	55.5	47.5	68.5	69.3	65.6	61.6	68.5	64.3
AA5	6.877	38.9	38.9	38.8	37.2	50.3	50.4	50.5	48.3	50.0	46.2
AA6	6.740	42.7	43.0	43.1	38.0	57.6	57.8	58.1	49.3	56.2	52.2
AB1	4.837	−80.2	−75.5	−75.5	−52.9	−95.2	−91.1	−90.7	−68.6	−96.8	−101.1
AB2	5.828	−65.3	−64.6	−64.2	−43.9	−76.6	−75.9	−75.6	−57.0	−79.7	−84.5
AB3	3.309	−60.1	−58.3	−67.4	−77.4	−79.1	−76.8	−90.8	−100.4	−82.8	−85.3
AB4	3.964	−86.7	−74.7	−77.0	−64.6	−107.7	−100.3	−102.0	−83.8	−106.2	−108.9
AB5	4.890	−52.1	−51.9	−51.6	−52.4	−70.0	−69.7	−69.7	−67.9	−69.7	−74.9
AB6	4.513	−52.8	−52.5	−52.2	−56.7	−70.4	−70.1	−69.4	−73.6	−70.2	−75.3
BB1	3.991	63.8	63.8	63.8	64.2	83.1	83.2	83.2	83.2	80.9	75.4
Adeninium chloride hemihydrate (ACH)
AA1	6.076	14.4	26.2	30.4	34.4	37.0	44.9	50.9	54.6	29.5	33.2
AA2	6.909	24.5	25.3	25.2	30.3	40.3	40.9	42.0	48.1	41.4	40.1
AA3	4.596	33.3	35.7	15.5	45.5	52.6	55.2	42.6	72.3	55.2	50.7
AA4	4.819	35.1	36.9	29.6	43.4	55.8	57.6	54.0	68.9	55.7	51.4
AA5	4.134	32.8	36.2	71.3	50.6	57.6	60.8	83.0	80.3	56.3	51.0
AB1	5.666	−63.1	−60.2	−58.0	−38.7	−88.1	−84.9	−82.8	−58.6	−86.2	−91.2
AB2	5.927	−45.3	−44.1	−43.1	−37.0	−58.7	−57.1	−56.3	−56.0	−61.0	−65.9
AB3	4.917	−54.3	−47.5	−47.5	−44.6	−82.5	−77.4	−75.5	−67.5	−81.1	−84.7
AB4	5.749	−49.8	−49.0	−48.5	−38.1	−69.0	−67.6	−67.7	−57.8	−72.7	−78.8
AB5	4.787	−50.7	−50.4	−48.9	−45.8	−76.5	−76.1	−75.0	−69.4	−76.7	−82.4
AB6	5.561	−40.4	−40.3	−41.1	−39.4	−57.5	−57.4	−57.9	−59.7	−57.7	−62.3
AB7	5.328	−42.2	−42.1	−43.0	−41.2	−60.5	−60.3	−61.0	−62.3	−59.4	−63.7
AW1	4.752	−7.7	−0.6	−1.2	4.2	−16.9	−12.2	−10.7	0.0	−14.8	−11.0
BW1	3.031	−23.8	−18.5	−18.6	−6.9	−19.1	−15.3	−15.2	0.0	−18.5	−14.3
BW2	3.659	0.6	0.7	0.7	−5.7	8.1	8.2	8.2	0.0	5.7	2.4
BB1	4.262	54.2	54.2	54.2	53.9	78.0	77.9	77.9	77.9	76.6	70.4
Guaninium Dichloride (GDC)
AA1	7.637	101.2	101.7	102.6	115.0	159.0	160.8	161.8	173.9	162.1	154.0
AA2	6.309	117.3	117.6	111.6	139.3	187.3	187.7	180.6	210.5	192.5	180.7
AA3	5.781	165.4	165.5	164.3	152.0	245.0	245.2	242.0	229.8	239.2	223.2
AB1	4.777	−123.1	−115.0	−115.4	−88.2	−184.3	−176.0	−176.2	−139.0	−181.4	−191.8
AB2	4.731	−98.5	−92.7	−92.7	−89.0	−153.3	−146.9	−147.1	−140.4	−153.6	−164.0
AB3	6.584	−80.3	−79.4	−78.6	−64.0	−119.0	−118.1	−117.6	−100.9	−119.5	−128.9
AB4	4.702	−123.0	−118.1	−116.0	−97.3	−177.9	−171.1	−170.4	−141.2	−169.9	−188.4
AB5	4.710	−95.4	−87.5	−87.2	−97.1	−144.6	−139.5	−138.5	−141.0	−151.5	−162.8
AB6	6.062	−95.3	−92.6	−91.1	−75.5	−136.3	−133.7	−130.8	−109.6	−138.5	−148.6
AB7	4.435	−91.5	−90.6	−88.3	−94.9	−144.9	−143.2	−138.0	−149.7	−144.2	−151.8
AB8	3.306	−102.1	−100.5	−111.9	−137.5	−151.5	−149.1	−160.4	−199.7	−156.0	−162.5
BB1	3.692	57.9	59.4	59.4	59.4	89.6	89.6	89.9	89.9	85.6	80.1

. As already discussed in [30], dimer interaction energies ($E_{tot}$) are in general much larger (in absolute values) than for analogous neutral dimers [59], and their sign is dominated by the net charges of the molecules, i.e., energies are negative for dimers of molecules having opposite signs of charge (AB type of dimers) and they are positive for dimers of molecules whose charges have the same sign (AA or BB type of dimers). Moreover, according to the experimental charge densities of the CC, ACH and GDC crystals, charges of nucleobase and chloride ions are not equal to the formal ones. In all the structures, chloride anions bear less negative charge, between −0.78(4) e and −0.88(5) e and, accordingly, nucleobase ions are less positively charged. Nevertheless, one has to assume that cations and anions extracted from crystal structures bear formal charges in order to apply methods such as UBDB, DFT–SAPT and PIXEL energy calculations which rely on monomer and dimer gas–phase (not periodic) approximation. Discrepancies in assigned molecular charges obviously have an impact on electrostatic energies computed for particular dimers, the $E_{es}$energies from experimental charge densities are significantly lower in absolute value, Table 2. However, as will be shown further on, these discrepancies do not influence the general conclusions which result from dimer electrostatic interaction energy analysis presented below.

Huge positive interaction energies for dimers containing nucleobase(s), meaning highly repulsive electrostatic interactions, seem to be intuitively in contradiction with supposedly attractive interactions identified on the basis of atom–atom contacts and topological analysis. Indeed, as stated in the introduction, it is commonly assumed that close atom–atom contacts between molecules bearing the same charge are forced to be present by attractive cation–anion interactions with neighbouring molecules [25]. However, if we take a closer look at the electrostatics, some understanding of the role of atom–atom close contacts in interactions between ionic molecules can be found. First of all, $E_{es}$ can be contrasted with energies computed applying Coulomb’s law to two point molecular charges placed at the distance at which centres of molecular masses are, see $E_{Coul}^{}$ in Table 2. Although point molecular charge is a very crude approximation of true charge density, it adequately explains the overall trends observed for $E_{es}$ values in relation to molecular charges and to distance. The trends are already visible for the set of the selected dimers presented in Table 2, but more appealing is Figure 2 in which energies for hundreds of dimers are plotted. The trends are qualitatively the same, regardless which method (PIXEL, DFT–SAPT, UBDB or even experimental or theoretical periodic data) was used to compute $E_{es}$ and $E_{Coul}^{}$ energies. The shorter the distance, the more dimers (but not all) exhibit energies differing from molecular point charge approximation. This is rather obvious since charge densities of nucleobases are far from being spherical (and point–like) and the closer the distance the more important the actual shape of the charge density of interacting molecules is. Interestingly, after subtracting the electrostatic interaction energy of molecular point charges ($E_{Coul}^{}$) from exact $E_{es}$, the remaining differences are values typical for dimers of neutral nucleobases in their crystals [59,63]. For some of dimers, the origin of the difference between the exact $E_{es}$ and the approximated one ($E_{Coul}^{}$) can be explained simply by the fact that higher electric multipole moments of nucleobase molecule(s) (higher than point charge) have important contributions. The inclusion of higher moments in energy calculations almost eliminates the difference, see $E_{mtp}^{mol}$ in Table 2. These are, for example, dimers AA5, AA6 and AB6 in CC and AA2, AB6, AB7 and BW2 in ACH. For the rest of the dimers, charge distribution of nucleobase(s) has to be described at atomic (not molecular) levels of detail to obtain energies close to the exact $E_{es}$ values. Nonetheless, for some of the dimers multipole moment approximation of atomic densities is still sufficient. These are AA4, AB2 in CC, AB4, AB5 in ACH, and AA1, AA2, AA3, AB3, AB7 in GDC, see $E_{mtp}^{atom}$ in Table 2. In the remaining dimers, the overlap of molecular (atomic) densities is large enough that the charge penetration effect on electrostatic energies cannot be omitted and exact integration over continuous charge density is necessary to obtain proper values of electrostatic interaction energy.

In many dimers, the molecules are oriented one to another to enhance electrostatic attraction, i.e., the exact $E_{es}$ values are less positive (for AA type of dimers) or more negative (for AB type of dimers) than $E_{mtp}^q$ energy resulting from molecular point charges. Large enhancement is, for example, observed for dimers AA1, AA2, AB1, AB2 and AB4 in CC, AA1, AA2, AA3, AA4, AA5, AB1, AB2, AB3, AB4, AB5, AW1, BW1 in ACH, and AA1, AA2, AB1, AB2, AB3, AB4 and AB6 in GDC. The largest enhancement due to strong hydrogen bonding was noted for the AA1 dimers in CC and ACH and described in the details in the first paper of the series, [30]. The other dimers usually contain strong hydrogen bonds as well or π...π interactions (as judged from distances and values of ρ(r) at BCPs) or several moderate ones, but this is not a rule. There are dimers for which, despite the presence of strong hydrogen bonds, the overall $E_{es}$ does not differ much from formal molecular charge approximation, so long range electrostatic interactions between the remaining parts of molecules seems to be less favorable. This is, for example, dimer AB5 in GDC. In some dimers, the difference even points towards a more repulsive type of electrostatic interactions than molecular point charge would suggest. These are usually dimers in which the π… π or π…Cl⁻ type of contacts were identified. Strong repulsive contribution is, for example, present in dimers AA4, AB3 in CC, AA3, AB8 in GDC, and BW2 in ACH. So molecules in these dimers are not oriented to each other in the optimal way from a molecular charge point of view, or the molecular charge is a far too simplified representation of charge densities for dimers having so many close atom–atom contacts.

3.3. A Missing Link between Topology and Energy Analyses

In the view of the above analysis, it is hard to find a consensus between the single atom–atom interaction and the whole molecule interaction views on the ionic crystal structures in this study. It is not true, for example, that the presence of strong hydrogen bond is always accompanied by lower repulsive electrostatic interactions between same charge molecules. It is also very hard to find any relation between exact $E_{es}$ (or $E_{tot}$) energies and $E_{EML}$ energies, even if analysis is performed on subgroups of like dimers (Figure 3a and Figure S5 in SI)). The lack of strong correlation between $E_{tot}$ and $E_{EML}$ energies we already observed in our previous studies [59] and was proven for large set of dimers by Spackman [64].

One more very interesting observation can be made regarding charge penetration. Whenever charge densities of molecules overlap significantly, electrostatic interactions between them cannot be quantified from point multipole approximation and integration over charge densities is necessary to obtain the exact values of $E_{es}$. The difference between exact values and values from point multipole approximation is often called penetration energy ($E_{pen}$). Although all the energy differences discussed above ($E_{es}$ and any $E_{mtp}$ or $E_{Coul}$) should, in fact, be called penetration energy, we choose to focus only on the difference between exact $E_{es}$ energy and energy resulting from point atomic multipole approximation ($E_{mtp}^{atom}$). We found such a defined $E_{pen}$ to be a very useful quantifier of the energetic consequences of electron density overlap. Penetration energy is usually negative, therefore, it enhances electrostatic attraction. It is present for short contacts, usually for those at distances around the sum of van der Waals radii and closer [65]. It appears that for all analysed structures, the $E_{pen}$ of a dimer correlates quite well with the sum of $E_{EML}$ for all BCPs encountered between molecules from that dimer; Figure 3b.

The correlation of $E_{EML}$ with $E_{pen}$ was found for all the three: experimentally derived crystal densities, crystal densities simulated by periodic quantum mechanics calculations and for crystal densities approximated by a sum of isolated molecule densities built from UBDB. However, it is stronger and the proportionality factor is closer to one for periodic models (experimental or theoretical) (Figure S6 in SI). Periodic models takes into account charge transfer between molecules and molecular charge density polarization occurring due to intermolecular interactions in the crystal as discussed in details in [30]. The UBDB model does not take into account these phenomena and can be named a procrystal charge density, i.e., a model of crystal charge density being a simple sum of charge densities of isolated molecules. This is similar to a definition of a promolecule being a simple sum of charge densities of non–interacting atoms. The UBDB model leads to closer to zero values of $\rho {\left(r\right)}_{BCP}$ and of $E_{EML}$ (Figure S4 in SI), as well as of $E_{pen}$ (Figure S7 in SI) when compared to experimental data. Apparently, $E_{pen}$ from UBDB differs more from experimental values than $E_{EML}$ leading to much larger than one $E_{EML}$/$E_{pen}$ ratio.

It is also necessary to investigate the influence of experimental errors and simplifications resulting from Abramov expression [62] and multipolar modelling on the observed relation. It is clear, for example, that Abramov expression (allowing to estimate values of the virial density $V{\left(r\right)}_{}$ at bond critical points of intermolecular interactions only from values of $\rho {\left(r\right)}_{BCP}$ and ∇²ρ(r)_BCP) may lead to errors of about 2–3 kcal mol⁻¹ for strong interactions, when compared to the values of $V{\left(r\right)}_{BCP}$ computed exactly from the wave function, see the E_EML plot in Figure S4 in SI.

Another interesting relations may also be found. For example, $E_{pen}$ correlates with $\rho {\left(r\right)}_{BCP}$ alone analogously to the $E_{pen}$ and E_EML, though not as strong. It appears that properties of BCPs, single point descriptors of interactions, represent very well the overall enhancement of electrostatic attraction due to the continuous nature of electron distributions of atoms and molecules [66,67].

4. Conclusions

The present work, following the work in [30], further deepens our understanding of intermolecular interactions among protonated nucleobases and chloride counter ions in the solid state. Moreover it gives some general remarks on the role of atom–atom contacts and charge penetration in understanding how charged (and neutral) molecules interact with each other.

We show how the overall crystal architecture of the studied compounds is governed by the tight packing principle and strong electrostatic attractions between molecules of opposite charges and strong repulsions between molecules with the same charge. Ultimately, minimization of the energy of the whole crystal decides the distance at which ions are placed and how they are oriented in relation to each other. Cations in the studied structures are larger and disc-like when compared to almost spherical anions. Small chloride anions generate stronger electrostatic field than bigger, single-protonated nucleobase cations having the same absolute charge. Thus, chlorides force nucleobases to form bidentate hydrogen bonded base pairs or π…π stacked dimers. In the case of adeninium, which is bigger than cytosinium, cations cannot be squeezed as much as cytosinium and neutral water molecules are incorporated into the crystal. Water molecules fill the space which cannot be filled by small chloride anions without disturbing the one-to-one balance between cations and anions. For doubly protonated guaninium, the relationship between the size and net charge of one guaninium dication, and the size and net charge of two chloride anions is well balanced, thus there is no need for guaninium to form bidentate base pairs and no need to include any additional neutral molecule into the crystal to dilute the charge. Hence, the studied crystals illustrate the idea of “Coulombic compression” as introduced by Braga et al. [25].

Interactions of molecular monopole moments are the major contributors to the $E_{es}^{}$ energies. However, the studied ions have shapes which are far from point–like and, in the case of nucleobases, far from being spherical. In addition, charge is not uniformly distributed inside an ion. Higher multipolar moments of molecules and of particular atoms, as well as the overlap of molecular charge densities matters. Thus, $E_{es}^{}$ energies for many ionic dimers in the studied crystals differ from the ones predicted solely from net molecular charge. The difference is often larger than 20%, and is sometimes even over 50%, and in one extreme case, i.e., dimer AA1 in CC (the trans sugar/sugar edge cytosinium base pair), exact $E_{es}^{}$ energy has a negative sign, indicating attractive electrostatic interactions between two nucleobase cations. In addition to this extremum, many other ions in direct contact are oriented in relation to each other in a way which strengthens their attractive electrostatic interaction or weakens strong repulsion, but not all of them. Usually, short hydrogen bonds or close stacking contacts are the source of electrostatic attraction contributions, but this is not a rule for the structures included in this study. Thus, these crystals exemplify the idea of “electrostatic attraction region” (EAR) as defined by Mata and co-workers [13], but also show that only certain dimers strictly follow this definition which evokes the overall minimization of cation–cation (or anion–anion) repulsion by EAR.

Despite their overall repulsive interactions, we observed numerous bond critical points (BCP) and bond paths between many nucleobase cations in the studied crystals using topological analysis of crystal electron densities (experimental or theoretical). Properties of BCPs do not differ much from those of neutral base pairs [63] as was observed for other systems [13] as well. Moreover, interaction energies computed from electron densities and their Laplacians at BCPs according to the Espinosa–Molins–Lecomte approach (E_EML) are negative, suggesting attractive interactions between cations. The values of E_EML are of similar magnitude for all short hydrogen bonds regardless if they are between nucleobase cations, nucleobase cations and water molecules or nucleobase cations and chloride anions. However, taking into account all the different types of interactions present in studied crystals, we have not observed any general relation between $E_{EML}^{}$ energies (their sum for particular dimers, to be exact) and $E_{es}^{}$ energies of dimers, similarly as was pointed out by Spackman [64] as well. On the other hand, we found strong correlations between $E_{EML}^{}$ energies and charge penetration contributions to the $E_{es}^{}$. In our opinion this is very significant observation which has never been reported yet (to our knowledge) in literature, although some discussion somehow envisioning the existence of this kind of relation is present in [68]. Our observation shows the quantitative link between the atom–atom close contact approach to structural analysis and analysis based on interactions between whole molecular charge densities. Atom–atom interactions identified using van der Waals sphere overlap concept, or better, topology analysis of crystal electron density, can be viewed as a manifestation of local electrostatic attraction of electron density of one atom with the deshielded nuclei of the other. Atom–atom contacts, especially in the context of ionic crystals, pinpoints the details of how local fragments of molecular charge densities interact with each other and what are the channels which strengthen electrostatic attraction through the short range effects related with electron density overlap.

We may speculate further that the observed relation between $E_{EML}^{}$ and $E_{pen}^{}$ has its roots not in classical electrostatic interactions but in other phenomena related to electron density overlap, like electron exchange. Pendas et al. [69] showed, for example, that the bond paths correspond to privileged channels of electron exchange and are a sign of a local lowering of exchange energy. Does $E_{EML}^{}$ quantify this change in energy?

Indeed, BCPs are important single point descriptors of interactions between whole molecules because they bear quantitative information regarding the overlap of molecular charge densities and its contribution to interaction energy. The charge density based approach for computing intermolecular interaction energies and the atom–atom approach to analyse non-covalent interactions do complement each other, even in ionic systems. Only by using these two approaches combined with structural analysis a comprehensive understanding of how molecules interact with each other in a supramolecular assembly, like a crystal or a macromolecule, can be achieved.

Analyses emphasising and quantifying the role of charge overlap might especially be important for the systems under high pressure, were chemical intuition related to the meaning of atom–atom contacts and their relations with interaction energies no longer works properly. For example, there are crystals in which not all atom–atom contacts shorten with rising pressure, some get longer [28,70].

Nevertheless, further investigations are necessary to fully understand physics behind the observed phenomenological relation between $E_{EML}^{}$ computed at BCPs and $E_{pen}$ for entire dimer.