Energy, orbital and structural stacking landscape of a purine homodimer system

The multidimensional study, combining the extensive calculations of potential energy surfaces for the parallel-displaced configurations and methods such as energy decomposition and natural bond orbital analysis, has been carried out. The resulted data give an energy, orbital and structural landscapes of this biologically essential system. The balance of the two energy sources, electrostatic and dispersion, is clearly visible. The obtained results, taken as a whole, provide an insight into the hierarchy of intermolecular interactions in the purine system, together with their sources.


Introduction
Noncovalent interactions (NCI hereafter) involving aromatic rings are pivotal in many areas of science. They are engaged in a vast array of supramolecular architectures playing a critical role in many chemistry-related fields such as the structures of biomolecules, molecular recognition, the reactivity of molecules, nanoengineering and in crystal engineering concerning systems being a collection of a discrete number of chemical units [1][2][3][4][5][6][7][8][9][10]. All the above strongly refers to N-heterocyclic compounds. Such systems, because of the presence of a nitrogen atom in their ring frameworks, are more abundant with π-electrons. Due to this, they can participate in different interaction patterns (e.g., π-stacking and C-H···N hydrogen bonds), which is often a key to their function [11][12][13][14][15][16][17][18]. Hence, N-heterocycles are useful recognition elements in many biological systems, and they are becoming a crucial component in most known drug molecules [14,[19][20][21][22]. Though the study of stacking interaction has mostly focused on systems containing benzene rings [23][24][25][26][27][28][29], the number of such studies on heterocyclic ring systems has been quickly increasing [11,14,16,[30][31][32][33][34]. Crystal engineering is based on the knowledge of intermolecular interactions. Because of that, there is no doubt that for it to be fully understood, it is necessary to know intermolecular interactions to the highest possible extent. The knowledge can benefit scientists who project and build new crystal systems as well as to those who analyze them.
Since a pyridine dimer is a prototypical example of π-π stacking in N-heterocyclic systems (as benzene dimer is an example of π-π stacking in general), it has been a subject of rigorous theoretical studies [11,15,31,32,[35][36][37][38][39][40][41][42]. In most cases, the studies have focused on the determination of system energy as a function of selected geometrical parameters, though some other approaches, as full geometry optimizations starting from various unsymmetrical initial dimer arrangements, have also been applied [43]. What is more, to enable a better understanding of the role of the molecular geometrical dependencies on the intermolecular interaction of the stacked pyridine dimer, the accurate potential energy maps of it, based on several geometrical parameters, have been calculated [11]. The same has been done for a pyrrole dimer [14] and, with some limits, to other simple aromatic N-heterocycles (e.g., pyrazine, pyrimidine, pyridazine) with benzene [44]. Recently, NCI between homodimers of selected heterocycles (pyrrole, furan, thiophene, pyridine, pyridazine, pyrimidine Electronic supplementary material The online version of this article (https ://doi.org/10.1007/s0021 4-020-02668 -y) contains supplementary material, which is available to authorized users. and pyrazine) have been studied to explore their interaction energy surfaces employing the Möller-Plesset second-order perturbation theory, coupled with small Gaussian basis sets (6-31G* and 6-31G**) with specifically tuned polarization exponents [45]. By using this approach, many vital conclusions have been deduced. It has become evident that the exploration of NCI landscapes may provide a new insight into NCI in aromatic systems [11,15,45]. As for more complex N-heterocycles, like purines, such rigorous studies, to the author's knowledge, have not been performed. The most studied systems are those found in nucleic acids, i.e., adenine and guanine [46][47][48][49][50][51]. The performed studies have mostly focused on analyzing the stationary geometries. The study of potential energy surfaces has been limited to learning potential energy curves showing dependencies of the system energy on twist or displacement [49,52]. In some studies, however, the potential energy surfaces have been examined by the molecular dynamics/quenching technique [53,54]. Then, the most stable and populated structures were fully reoptimized at the correlated ab initio level [53,54]. Because of the above, as the focus of this study, the extensive analysis of potential energy surfaces (PES) of model purine-purine dimers, combined with orbital and structural analysis, has been selected. Accurate PES maps of a purine system would enhance the knowledge of the complete understanding of the stacking interaction of this biologically fundamental structural unit. Unlike analyzing simple potential energy curves, one could see a broader picture of what occurs and learn how the purine dimer stability changes along with various geometrical parameters. Accurate multidimensional PES maps, as opposed to the data based only on a small group of systems, could find a vast range of applications in chemistry. Such maps would appear to be a guide for qualitative predictions of binding interactions. The obtained data may be related to other aromatic systems, especially those purine-based like adenine and guanine. Hence, new aspects of base stacking in DNA may be discovered. The results may also help to learn more of the role of the position of nitrogen atoms and C-H/N-H bonds (linked to the mutual position of the interactions molecules) on the system energy as well as determine the geometrical boundaries in which the interaction energy stays the same or changes negligibly. For the above reasons, the performed research may translate into a better harnessing the full potential of NCI (e.g., in predicting the structures of supramolecular systems and in molecularly targeted therapies). By combining the data from Cambridge Structural Database (CSD) [55] with the extensive energy analysis, and the use of methods such as energy decomposition and natural bond orbital (NBO) analysis, the performed study can disclose much deeper crucial aspects associated with the physics of NCI (e.g., the actual importance of π orbitals into intermolecular interactions in stacked systems).
The accurate calculations of NCI energies require an appropriate description of dispersion forces. This can be achieved by using advanced wave function theory (WFT) methods. Unfortunately, obtaining PES maps in this way remains a challenge. This stays true even for small aromatic systems since the advanced WFT methods are very expensive and cannot be applied for routine calculations [56]. However, an alternative is the use of density functional theory (DFT) methods whose cost is significantly lower. Over the recent years, many new density functionals designed for the energy calculations of non-covalently bonded systems have been developed [57][58][59][60]. Dispersion forces may be included by adding a dispersion correction term. This method is straightforward and requires no additional computing resources. Some dispersion-corrected density functionals with D3 or D3(BJ) correction lead to the results that are close to those from advanced WFT methods for dispersion-dominated NCI [57,58,61].

Methodology of theoretical calculations
Model purine homodimer systems were constructed in the way shown in Fig. 1. In each of them, the purine frameworks were parallel to each other. The energy calculations were carried out as a function of four geometrical parameters: the distance between the centroids of the purine frameworks (d), the twist angle (α) between the purine frameworks, the angle between the line connecting the purine framework centroids and the normal line to the purine framework in which the connecting centroid line starts (β) and the angle associated with the rotation of one purine monomer around the other (γ). The calculations were performed employing Gauss-ian09 Rev. D.01 [62]. The starting monomer was a purine molecule with geometry optimized at MP2/aug-cc-pVDZ level, where MP2 stands for second-order Møller-Plesset perturbation theory [63]. The energies of all dimers were corrected for basis set superposition error using the counterpoise method [64]. Interaction energies were calculated by subtracting the energy of the two monomers: For γ = 0°, the shift of one ring along the plane of the other was following the direction of the vector formed between the other framework centroid and the carbon atom labeled as 5 (Fig. 1, the vector beginning was the position of the framework centroid).
An appropriate DFT functional was found by performing the energy scan of the selected purine-purine homodimers ( Fig. 2). The scan was based on CCSD(T)/CBS level of theory using aug-cc-pVDZ and aug-cc-pVTZ basis sets [65] to obtain the MP2/CBS limit energy through the extrapolation scheme as described by Helgaker et al. [66]: where X = 2 and Y = 3 (for the used aug-cc-pVDZ → aug-cc-pVTZ extrapolation). E (X) and E (Y) are the MP2 energy values calculated using aug-cc-pVDZ and aug-cc-pVTZ basis sets, respectively, for E (X) and E (Y) . The difference between CCSD(T) and MP2 energy has insignificant (unless a too small basis set is used) basis set dependence [67]. Therefore, to result in CCSD(T)/CBS energy, the difference between CCSD(T)/aug-cc-pVDZ and MP2/aug-cc-pVDZ was added to the extrapolated MP2/CBS energy. The purine-purine energies were calculated as a function of d. The distance d was varied from 3 to 9 Å with 0.5 Å increments. The angles β was set as 20° and α and γ were set as 0°. The value 20° for β angle was based on the author's experience. It is also generally accepted for π···π bonding interactions. [11,14,15] The calculations were then repeated using two density functionals: B3LYP-D3(BJ) and M06-2X. B3LYP-D3(BJ) is a dispersion-corrected DFT with D3 version of Grimme's dispersion with Becke-Johnson damping [58]. It also successfully described NCI in a system previously studied by the author [11]. M06-2X is a popular functional which performs well in the calculation of NCI between aromatic molecules [59,60]. Both functionals employed cc-pVDZ basis Interaction energy curves created on the basis of data obtained for the selected purine homodimers. The curves were created using splicing interpolation. CBS stands for a complete basis set extrapolation set. All CCSD(T) calculations have been performed using Psi4 Package [68], and they were utilizing density fitting (DF) algorithm [69]. The best functional happened to be B3LYP with D3(BJ) correction (Fig. 2), and it was used to perform all the further calculations. The purine homodimer energies were calculated as a function of four selected parameters (Fig. 1). The parameter d was varied from 3 to 9 Å with 0.1 Å increments, α was varied from 0 to 345° with 15° increments, β was varied from 0° to 90° with 10° increments, and γ was changed from 0 to 180° with 45° increments. Calculations for γ exceeding 180° were not needed since the obtained geometries were covered in those created for smaller γ. For γ = 0°, the lateral shifting (related to the increase in d at the given β) of one ring along the plane of the other was following the direction of the vector created between the other framework centroid and the carbon atom labeled as 5 (Fig. 1, the vector beginning was the position of the framework centroid).
Based on the gathered data, the contour plots representing energy as a function of α and β and as a function of d and β (Fig. 3) were created using OriginPro 2020 [70]. Some intermolecular interactions in the investigated purine-purine systems (configurations corresponding to the found energy minima and those related to PES map areas for which the changes in energy are negligible) were analyzed utilizing localized molecular orbital energy decomposition analysis (LMO-EDA) [71]. It is a method, which is an extension and modification of the techniques developed by Kitaura and Morokuma, Ziegler and Rauk, and Hayes and Stone [72][73][74]. LMO-EDA is easily used with both DFT and WFT approaches. It has been successful in the study of the intermolecular interactions (e.g., hydrogen bonds, π···π and C-H···π contacts) in many The PES maps are shown together with the corresponding values of rotation angle between the monomers (a) and the corresponding distances between purine frameworks centroids (b) systems such as pyrrole and pyridine homodimers [11,14], indole-cation-anion complexes [75], aromatic units of amino acids with guanidinium cation [76], hydrogenbonded complexes of serotonin [77], radical conjugated systems [78] and C-H/π complexes in water [79]. This method decomposes the interaction energy to electrostatic (E ele ), exchange (E ex ), repulsion (E rep ) and polarization (E pol ) terms, and they are calculated based on single-determinant Hartree-Fock (HF) wavefunctions. The correlation term (E corr ) is derived from a supermolecule approach by the use of a correlation method (e.g., MP2 and CC), and it equals the difference between the correlation method energy and HF energy [71]. In the performed study, the contribution of correlation energy was estimated from B3LYP-D3BJ/cc-pVDZ calculations. It was calculated as a difference between B3LYP-D3BJ/cc-pVDZ energy and HF/cc-pVDZ energy. LMO-EDA calculations were carried out with GAMESS [80,81] software employing HF/ cc-pVDZ method. As implemented in GAMESS, in LMO-EDA the counterpoise correction for basis set superposition error is used.
The study of potential orbital interactions between the purine units covered the model purine homodimers, which configurations corresponded to the found energy minima, both in a function of α and β. In this study, NBO analysis was performed at M06-2X/aug-cc-pVTZ level. Homodimer geometries corresponding to local energy minima (in function of β) were also analyzed, at M06-2X/aug-cc-pVTZ level of theory, in terms of Hu, Lu, and Yang (HLY) charge-fitting method [82]. According to the literature, it gives better results than those obtained from commonly used CHelpG scheme [82,83]. In the cases mentioned above, M06-2X density functional was selected due to its excellent performance for application in the general chemistry of main group elements [59,60].
For the statistical analysis, the compounds included in the Cambridge Structural Database (version 5.41, November 2019) containing purine frameworks being close to parallelism (the maximum dihedral angle between these ring planes was 10°) and with the distances between the ring centroids being in the range of 3-9 Å were selected. Based on the obtained data, the pairs of following parameters were determined: the purine framework centroid distance (d) and the angle between the line connecting the purine framework centroids and the normal line to the purine framework in which the connecting centroid line starts (β) (Fig. 1). The obtained pairs were then analyzed using twoDBinning tool as implemented in OriginPro 2020. The bin sizes were 0.5 and 10, respectively, for d and β. The occurrence frequency obtained in this procedure data pairs (d, β) was compared with the obtained PES maps.

Results and discussion
Out of the two selected density functionals, the best one happened to be B3LYP with D3(BJ) correction. It gave superior results comparing to M06-2X (Fig. 2). For it, the mean deviation from CCSD(T)/CBS results was equal to 0.31 kcal/mol. However, the worst outcome for it was associated only with a sandwich configuration (with no lateral shifting). Excluding this configuration, the mean deviation is equal to only 0.10 kcal/mol. For M06-2X, the mean deviation is much higher, and it is equal to 1.17 kcal/ mol. Additionally, the differences are also seen concerning the binding regions. M06-2X predicts well the geometry associated with the minimum energy, but unfortunately, it gives non-binding interaction energies for dimers for which CCSD(T)/CBS gives binding interaction energies. Because of this relatively significant underestimation in binding energy (Fig. 2), using M06-2X might lead to the formulation of inappropriate conclusions on the studied system. B3LYP with D3(BJ) correction reflects well the binding regions predicted by CCSD(T)/CBS. For the studied system, it was also significantly faster, almost two times, than M06-2X. The computational time ratio for the selected functionals is 1:1.8, where the values correspond to B3LPY-D3(BJ) and M06-2X, respectively. Taking it all into account, B3LYP with D3(BJ) appears to be the most suitable method for the analysis of interaction energy in the purine homodimer system.
Multidimensional energy scan allowed to create an energy landscape of the purine system. One can quickly see how the purine dimer stability changes along with various geometrical parameters. Based on the obtained maps, the optimal configurations of the purine homodimer may also be identified (Fig. 3). One can also find that the binding boundaries (related to the areas with negative interaction energies values) are quite broad. They correspond to a wide range of selected geometrical parameters. Even the minimum d, equals to 3 Å, is enough to provide binding interaction energies. In the beginning, with the increase in d, the interaction energy between the monomers decreases. This encourages the purine molecules to bind more strongly to each other (Fig. 3). Yet, with the further increase in d, the interaction energy starts to get larger. Generally, the most stable configurations can be found in the "valley" extending from d in the range 3 and 4 Å for β = 0° and changing gradually into the values in the range 5 and 9 Å for β = 90°. An optimal system configuration cannot be reached only by the separation of monomers. Other parameters need to be changed as well. Because of this, a sandwich configuration, though energetically stable as it may seem (Fig. 3), is not associated with an energy minimum. For this configuration, the electrostatic energy term, though negative, is higher compared to others (Fig. 4). To decrease it-make it more optimal-one of the monomers needs to be shifted along the plane of the other. This shift decreases not only the electrostatic energy terms but also polarization and correlation energy terms, which, for β = 20°, are in their local minima (Fig. 4). One must keep in mind, however, that these local minima for the described energy term factors are set in the space of the configurations with their minima for the given β, d pairs and are a function of β (Fig. 3). The particular energy terms are discussed further in the text. The majority of the configurations with a minimum binding energy for the given d and β are found for γ = 0° which is the angle determining the rotation of one monomer around another (Fig. 3). Around half of the "valley's bottom" is associated with this γ value. The most energetically favorable configuration, though, is related to γ = 90°. For this geometry, the interaction energy between monomers is equal to -13.89 kcal/mol (Table 1) and this low value is associated with the presence of strong dimeric N-H···N hydrogen bonds (Fig. 5). For β = 90°, for which the interacting monomers lie in the same plane, the range of binding boundaries for d is relatively wide and it expands up to the calculation range (Fig. 3). The particular energy landscapes of the interaction energy values seen as a function of d, β pair are similar when they are observed separately for the respective γ. However, this similarity is somewhat associated with the boundaries of binding energies seen in terms of d and β. The interaction energy values for the respective γ differ sharply. Those differences, though, are the most significant for high β for which it is possible to engage the interacting monomers into the formation of N-H···N hydrogen bonds. In the case of such configuration, for β equal to 90°, both electrostatic and polarization energy terms are the lowest (Fig. 4). When it comes to "typical" stacking (for smaller β and d), the interaction values are similar for each γ, and they equal to around -7 kcal/mol (Supplementary Materials, Figure S1).

Fig. 4
Decomposed energy terms for stacked purine homodimers whose configurations correspond to energy minima found for the systems with the given β. The curves were created using splicing interpolation An intriguing landscape appears when the system energy is perceived as a function of α and β parameters (Fig. 3b). It is straightforward to see that there are several "deep" energy "valleys" for high β values (in the range of 70-90°). For much smaller β, in the range of 0 to 40°, there is a great energy "valley." It is not as deep as the "valleys" for higher β values, but its area is the largest. It is associated with the distance in the range of around 3-5 Å and covers stacked geometries. The interaction energy between monomers in those is similar, and it is equal to about -7 kcal/mol. The correlation energy term is slightly dependent on α (Fig. 6), and this is related to the interaction surface of the monomers-keeping Decomposed energy terms for stacked purine homodimers whose configurations correspond to energy minima found for dimers corresponding to configurations with the given α values and with β = 20°. The curves were created using splicing interpolation all other parameters the same and changing only α the interaction surface of the monomers is different when one monomer is twisted around its vertical axis. For the stacked geometries and β = 20°, those changes in correlation energy terms are in the order of around 2 kcal/mole (Fig. 6). Other terms are more dependent on α, especially for electrostatic energy term (Fig. 6). This is associated with the presence of nitrogen atoms in the purine frameworks. Their presence influences the specific charge distribution in the monomer. Hence, some dimer configurations are electrostatically much less favorable than others. The "valleys" for higher β are associated with the formation of hydrogen bonds. The creation of this type of bonding is feasible for several α. However, the most energetically stable configurations are observed for α equaled around 180°-220°, where strong dimeric N-H···N hydrogen bonds may be formed (Fig. 3,  Fig. 5). For those α values, the geometries associated with their energy minima correspond to such systems. In their case, both electrostatic and polarization energy terms are the lowest and are much lower from those terms found for configurations correspond to energy minima for other α (Fig. 7). Analyzing the graph showing the changes in correlation energy terms (related to the structures with energy minima for the given α), one can see with ease for which alpha the hydrogen-bonded structures are present. For those, the correlation energy terms are the highest and are similar to each other (Fig. 7). Monomers in such systems are coplanar or are close to coplanarity. Hence, the interaction surface is significantly smaller than for stacked units. In this landscape, perceived as a function of α and β parameters, there are no energetically unfavored configurations. For each α, β pair it is possible to find a system with negative interaction energy. The smallest binding energies correspond to α equaled 150°. Considering the individual γ angles, each energy landscape is similar when it comes to the region when two purine molecules are in a typical stacked configuration (β, in the range of 0° to 40°). Considerable differences appear for high β where two purine units are close to coplanarity. In this configuration, the formation of N-H···N hydrogen bonds is possible, and this formation is dependent on γ that is an angle determining the rotation of one molecule around another. For the respective γ, the interaction energies for some configurations are not binding. Those configurations are the ones with the C-H bonds pointing to each other (Supplementary Materials, Figure S2).
With the use of NBO, the above-described energy landscapes lead to a "quantitative picture" of orbital interactions which reveal their essential role in the intermolecular bonding. In the case of dimers whose geometries correspond to energy minima (Fig. 5), the importance of the interactions which are related to the overlap of the molecular orbitals is shown (Table 2, Fig. 8). For the configuration Fig. 7 Decomposed energy terms for stacked purine homodimers whose configurations correspond to energy minima found for the systems with the given α. The curves were created using splicing interpolation associated with a stacked geometry (Fig. 5), the most critical orbital interactions are the interactions between π bonding and π antibonding orbitals that are present in C-C and C-N bonds of the interacting monomers (Table 2, Fig. 8). It is the same for other stacked geometries. Unless other geometrical parameters are not energetically unfavorable (Fig. 3), the π → π* interactions are dominative for β, in the range of 0 to 60° (Supplementary Materials,  Table S1-S2). Hence, π → π* interactions are significant even for stacked geometries for which one would instead expect the dominative role of C-H···π or N-H···π orbital interactions. The same was observed in the case of a pyridine dimer system [11]. Considering stacked geometries, the π → π* interactions are not the only important ones. The others are the interactions engaging a lone pair (lp) of nitrogen atoms in which the lp donates the electron density to π antibonding orbitals or some Rydberg orbitals (e.g., carbon Rydberg orbitals or hydrogen Rydberg orbitals). The interactions between σ orbitals of C-C and C-N bonds and Rydberg orbitals of C atoms have also been found (Table S1-S2). With the increase in β, when it exceeds 60°, the π → π* interactions vanish and some other appear. The number of different orbital interactions rises. Mainly, the interactions involving a lone pair of nitrogen atoms are observed (Table S1-S2). For instance, one may find lp → ry* H interactions for which a lone pair of a nitrogen atom donates to Rydberg orbital of a hydrogen atom. Nevertheless, the ones of the most important are lp → σ* N-H in which the lp of a nitrogen atom, which is not connected with a hydrogen atom, donates electron density to an antibonding σ orbital of an N-H bond (Fig. 8). For β =90°, the stabilization energies of such orbital interactions are very high, and for some configurations, they exceed 10 kcal/mol (Table S1-S2). The changes in charge distribution confirm the redistribution of electron density. For β =90° and for a configuration for which strong dimeric N-H···N hydrogen bonds can form, the differences in charge distribution of the interacting monomers are easily seen compared to a non-interacting monomer (Figs. 1, 5). The nitrogen atoms which are considered donors become less negative and the carbon atoms, labeled 1 and 14 (Fig. 5), which are present between nitrogen atoms connected with hydrogen atoms and nitrogen atoms which donates, become less positive.  For stacked systems, such significant differences in charge distribution are not seen (Figs. 1, 5). The search of the CSD for structures containing purine as a structural unit resulted in 4037 hits. The number of structures seen in a function of the geometrical parameters d and β is shown in Fig. 9. The obtained data, seen graphically as the function mentioned above, reflects the results from the performed calculations (Figs. 3, 9)-the number of the respective structures gets more significant as d and β increase. The observed relationship is analogous to the changes in the interaction energy between monomers. However, despite this general similarity, some differences are present. The first obtained population maximum seen at β equaled around 30° and at d equaled around 4 Å corresponds well to the energy minimum seen on the calculated PES map (Fig. 3). It is associated with the stacked purine units being in a parallel-displaced configuration. In those structures, the purine units are not incorporated into large molecular framework systems. In most cases, they are a part of adenine sub-units (Fig. 10). The second population maximum (Fig. 9), which can be easily distinguished, is observed at β in the range of 45-55° and d within the range of 4.5-5.5 Å. It is associated with more structures than the one mentioned earlier, and it is not associated with the energy minimum on the respective PES map (Fig. 3). In the case of it, one can find many configurations with purine units being a part of an adenine sub-unit or an adenine derivative subunit (Fig. 10). This population maximum is mostly associated with arrangements corresponding to N-H···π (engaging -NH 2 group of adenine sub-unit) or C-H···π interactions. The third population maximum, at β in the range of 80-90°  [85], c AHEWIT [86] and d in the rage of 7-8 Å, involves structures containing mainly adenine and adenine derivative sub-units. Compared to the energy minimum on the calculated PES map (Fig. 3), it is shifted by about 1 Å toward more considerable distances. The reason is the formation of N-H···N hydrogen bonds in which the -NH 2 group and a nitrogen atom of pyrimidine or imidazole rings are involved (Fig. 10).
Looking at the structural stacking landscape based on the data found from CSD, one can see in its lower right quadrant a "peninsula of structures" extending from d equaled around 7-7.5 Å for the smallest β angles to d equaled 9 Å for β in the range of 40-50° (Fig. 9). According to the calculations, the area associated with this "peninsula" is energetically binding. However, the PES obtained from the calculation does not have an energy minimum in that region (Figs. 3,9). Because of that, the presence of such a large number of structures in this area is quite intriguing. The configurations corresponding to that region cover the purine frameworks that are not coplanar and with quite large distances between the interacting purine units. Some geometries are even in a sandwich conformation (for β = 0°). All of this makes the presence of this "peninsula of structures" in the CSD-based landscape even something exotic. However, analyzing the structures related to that region that were found in the CSD, it turned out that the respective purine units do not interact directly with each other. Between the interacting units, there is an auxiliary system, e.g., another aromatic unit or polyiodide anions (Fig. 10). In the case of other studied systems (pyridine-pyridine, benzene-benzene and benzene-pyridine), this phenomenon was not observed in their CSD structural landscapes [15].

Conclusion
In terms of binding boundaries, the PES maps are similar to those found for other homodimer systems (i.e., benzene, pyrrole and pyridine) [11,14,15]. For a purine dimer, opposed to those earlier mentioned, the maximum binding energy is related to a configuration where two monomers are coplanar (a hydrogen-bonded configuration). This confirms the relatively high number of structures containing coplanar purine units (Fig. 9), where those units interact with each other via hydrogen bonds (Fig. 10).
The landscape related to the balance of the two energy sources, electrostatic and dispersion, is clearly visible. The contribution of correlation energy to the energy seen as the sum of electrostatic and correlation terms is equal to around 60-70% for stacked geometries (Fig. 4). It seems to be a little smaller than observed for benzene and pyridine [11,14], but it is significantly higher than that found for pyrrole equaled around 50% [14].
The role of orbital interactions has been found, and their hierarchy in this system has been known. In the case of stacked monomers, the interaction related to π-orbitals overlap may provide additional binding force. The same was found for pyridine dimers [11]. That orbital overlap is important even for configuration for which one would instead expect the formation of N-H or C-H···π interactions. It is worth mentioning that for pyrrole, significant orbital interactions have not been observed [14]. The above reflects the complexity of the issue.
The geometrical boundaries for which the interaction energies are binding are very broad (Fig. 3). As a result, one may have a great adjustability in the crystal engineering processes. What is essential, however, is to take the effect of such interactions into account for systems that one may analyze, especially when coplanar ones are considered. The number of structures containing coplanar purine units is relatively high (Fig. 9).
The conducted research shows the complexity of the interactions between purine monomers. Though the presented results shed more light on the stacking phenomenon, there is still yet much to learn and discover, especially in terms of the hierarchy of intermolecular interactions and their sources. There is still no unified description of the factors that determine this chemical phenomenon in general.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.