Fractal nature of benzene stacking interactions

Benzene and other aromatic groups, as planar groups with π\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi $$\end{document} electrons cloud, tend to form stacking interactions which have an important role in various chemical and biological processes. In order to have a better insight in the nature of these interactions, we have performed a fractal analysis on patterns of electron density and electrostatic potential for two benzenes in stacking interaction. The calculated fractal dimension follows the trend of the calculated interaction energy for the interplanar distances of 4.0 to 6.0 Å, which partially coincides with the strongest attractive stacking interactions. The fractal dimension vs. energy dependences were fitted with the logistic curve, and the fitting coefficient was 0.96 up to 1.00. For the benzene stacking interaction energy, with a range of conformations and distances between two benzenes, DFT calculations at the B3LYP+D3/aug-cc-pVDZ level were performed with the TURBOMOLE software. The fractal analysis for electron density and electrostatic potential has been done by python scripting.


Introduction
Planar organic molecules which posses π electrons tend to form stacked multimers. There are several models explaining the nature of π -π interactions [1][2][3][4][5][6][7][8], and our intention is to contribute with another possible one. We sensed that the interaction energy of these systems must be somehow related to fractality. So we have analyzed patterns of distribution of a spatial molecular property related to interaction, such as electron density and electrostatic potential.
A fractal is an object that has the property of self-similarity and scale invariance. Its Hausdorff-Besicovich dimension is strictly larger than the topological dimension [9]. And its Institute of Chemistry, Technology and Metallurgy -National Institute of the Republic of Serbia, University of Belgrade, Njegoševa 12, Belgrade 11000, Serbia dimension is fractal, i.e., not an integer. Fractal dimension describes complex patterns (natural ones) which are hardly defined in Euclidean space. Determining the fractal dimension (fractal analysis) is used for pattern characterisation in different fields [10]. Self-similarity involves a pattern which is composed of the same but smaller copies. Because of this feature, fractals are invariant on the magnification scale. Euclidean dimension is an integer: 0 for point, 1 for line, 2 for plane... In plane, a fractal curve has fractal dimension between 1 and 2. Theoretical maximum is 2 when the fractal curve becomes so complex that it fills the plane.
Fractal dimension is calculated like this [9]: where N s is the number of self-similar pieces, s is the size of the pieces (or the linear scaling factor of the pieces to the whole), and F D is the fractal dimension that we want to calculate. It can only be strictly computed for the deterministic fractals. Natural objects practically do not exhibit deterministic self-similarity, but fortunately, they exhibit some statistical self-similarity. Therefore, the box-counting method is used [11]. This method is based on the iterative counting of boxes that contain at least one piece of the frac-tal pattern, for example, at least one black pixel pertaining to a black curve on a white plane. The count is performed for every scale factor, or box size. After that, the FD can be estimated from the least squares linear fit of log(N s ) versus log(s). For grey scale images, the differential box-counting method has been proposed [12][13][14][15]. For a two-dimensional image, the grey intensity would be the third dimension, so the fractal dimension would be between 2 and 3 for 2D images. We will use the differential box-counting method because it suits our patterns with varying intensities of electron density and electrostatic potential distributions.

Results
We have calculated the interaction energy for two parallel stacked benzenes with offsets and determined the fractal dimension for every image of electron density map intersection for the three orientations, a, b, and c as in [16] (Fig. 1), for offsets r ranging from 0.0 to 6.0 Å with the step of 0.5 Å. Previous calculations have shown that the strongest interactions are not the sandwich ones but the ones with significant offsets [16][17][18]. The interplanar distance R is 4 Å which corresponds to the strongest interactions for zero offset calculated by the method in this work. For all the three orientations, the fractal dimension increases with energy following very similar curve, and the fractal dimension minimum corresponds to the energy minimum (1.5 Å) (Figs. 2, 3, and 4). Fractal dimension dependence on energy exhibits sigmoidal shape, so we have chosen the logistic function: and fitted it to the obtained curves, because this function includes the exp(-x) term which is interesting in the light of statistical thermodynamics approach, if the x variable is energy. The obtained agreement is very good, represented by the coefficient of determination of the logistic regression k 2 . The obtained FD(E) functions are given in the figures captions. The functions parameters are similar for the different orientations.
Furthermore, we have performed the same analysis on the electrostatic potential maps (Figs. 5, 6, and 7). For the electrostatic potential, the patterns were treated as grey images as well, but they are shown in the red-white-blue palette in order to represent regions of positive and negative potentials. Electrostatic potential as well shows a good agreement, with very similar k 2 values, but in contrast to the electron density, the fractal dimension minima are offset by 0.5 Å in relation to the energy minima for the orientations a and b, although the minima are very slightly pronounced. Electron density might be more straightforwardly related to fractal dimensions because it is a physical quantity that is more fundamental, not derived as electrostatic potential which is a function of electron density.
In order to elucidate more in details the nature of the stacking interactions, different energy terms such as electrostatic and dispersion energy were compared separately to the fractal dimension. Table 1 shows the k 2 values of logistic regression for the energy decomposition. All the energy terms exhibit as good fit as the total energy, k 2 from 0.98 to 1.0, except for the electrostatic energy.
As the benzene interaction energies are significant for the interplanar distances R other than 4.0 Å, we have also performed the fractal analysis for R values from 2.0 to 8.0 Å with the step of 0.5 Å and offsets r varying from 0.0 to 6.0 Å for each R. The fractal dimension does not always Fig. 1 The three orientations of the parallel stacked benzene dimers used for the calculations [16]: offset along the C-H bond (a), staggered offset (b), and offset perpendicular to the C-H bond (c). Offset r is the projection of the distance between the centres of the two benzene rings to the x-y plane, and the interplanar distance R is the z-axis projection of the distance between the centres of the two benzene rings. The electron density or electrostatic potential map intersection (represented in colour) is the x-y plane in the middle between the benzene ring planes; its z coordinate is R/2  Table 2). The best logistic regression is for R of 4.0 to 6.0 Å which partially corresponds to the strongest attractive stacking interactions. The calculated interaction energies for all the displacements R are represented in Table 3.

Conclusion
All the calculated electron density and electrostatic potential patterns for benzene stacking interaction are actual stochastic fractals as their calculated fractal dimensions have Unexpectedly, the fractal dimension increases with the energy, so the conformation with the strongest interaction exhibits the lowest fractal dimension. One would expect the highest fractal dimension, i.e., the most complex space filling for the most attractive interaction.
Energy decomposition analysis shows as good logistic fit against the fractal dimension for all the energy terms separately, as in the case of the total energy, except for the electrostatic energy.
Even though we have not considered the surroundings of the benzene molecules, or other aspects as entropy or time, this finding might still be a significant clue for deeper understanding of the π -π interactions.
The electron density and electrostatic potential maps were defined in a way that they encompass the whole system of the two molecules with padding of at least 5 Å, and the molecules L is the image length in pixels. The box sizes s are ranging from 2 to L. G is the total number of colours, G = colour max − colour min + 1. Colours are ranging from 1 to L and are converted to integers. For comparison, all the images are of the same size and resolution. The box counting N s is then plotted against box sizes, and the fractal dimension is determined as negative slope of the logarithmic plot log(N s ) versus log(s). The box size ranges from the size of the whole image and then splitting the image 2 times consecutively until the size of 1 pixel is reached. Therefore, we adjusted the image to be squared and with the size which is the exponent  probe resolution of 1024 x 1024 gave as good agreement of energy and fractal dimension. If a logarithmic plot is linear, then the pattern is a fractal. The error of fractal dimension calculation is 1-coefficient from linear fitting of the logarithmic plot. All our images have very good regression, i.e., small fractal dimension error, so they are indeed stochastic fractals.
All the scripts for fractal analysis, plotting, and curve fitting were written in python [31,32]. The VMD software was used for visualisation [33].