Theoretical study on prismatic (N6)n (n=16–35) molecules

Journal of Molecular Graphics and Modelling 96 (2020) 107508

Contents lists available at ScienceDirect

Journal of Molecular Graphics and Modelling journal homepage: www.elsevier.com/locate/JMGM

Theoretical study on prismatic (N6)n (n¼16e35) molecules Qian Guo a, Bing He b, Hongwei Zhou b, * a b

College of Mathematics, Chengdu Normal University, Chengdu, 611130, PR China College of Chemistry and Life Science, Institute of Functional Molecules, Chengdu Normal University, Chengdu, 611130, PR China

a r t i c l e i n f o

a b s t r a c t

Article history: Received 20 September 2019 Received in revised form 4 November 2019 Accepted 4 December 2019 Available online 14 December 2019

As a continuation of our study on the stability of pure nitrogen molecules, twenty members of the (N6)n (n ¼ 16e35) molecular sequence with D3h or D3d symmetry alternatively are studied in this work. The structures and energies are examined at the B3LYP/cc‒pvDZ computational level. “Natural Bond Orbital (NBO)” and “Atom In Molecule (AIM)” analyses are performed to investigate the bonding properties and the electronic topologies of the prismatic molecules. The van der Waals forces or the dispersion interactions are identified to be the dominant stabilizing factor for large prismatic nitrogen cages. The van der Waals forces interweave with covalent bonds to form a dense network which tightly bind nitrogen atoms at grid points, and consequently, to make the prismatic (N6)n molecules stable. The mechanism of generating van der Waals forces in the prismatic molecules is the focus of this paper. Geometrically, all the polygons on the prismatic surface are boat hexagons except the two caps. The van der Waals forces or dispersion interactions occur only between the nitrogen atoms at positions 10,40 of the boat hexagon on the surface. The atom‒atom proximity effect of nitrogen atoms at the positions 10,40 leads to the van der Waals interaction. Topologically, van der Waals forces are generated by partial overlap of sp3 hybrid orbitals of lone pairs of nitrogen atoms at positions 10,4’. These large prism‒shaped nitrogen cages may also be novel beeline nanotubes which is environmentally friendly. This work is expected to open up a bright prospect of the study and applications of nitrogen nanofibers. © 2019 Published by Elsevier Inc.

Keywords: Atom‒atom proximity Molecular sequence van der waals forces Prismatic Boat hexagon Chair hexagon unit

1. Introduction As a potential candidate material with high energy density, pure nitrogen compounds have been studied extensively by theoretical chemists from the 1990s to the early 2000s [1e5]. Carefully reading these studies, it has been found that their discussions mainly focus on the thermodynamic stability of these pure nitrogen compounds, since only a compound with sufficient thermodynamic stability can be synthesized, and as an efficient propellant and high explosive, its decomposition reaction Nx / (x/2)N2 is controllable. The early ab initial studies on thermodynamic stability of pure nitrogen molecules were restricted only among the relatively small molecules [6e13]. In order to conduct computational research on larger pure nitrogen molecules (such as N60), a semidempirical calculation which requires less computer resources was applied [14]. With the progress of computer technology, larger pure nitrogen molecules with various structural forms have been studied in a computational

* Corresponding author. E-mail addresses: [email protected] (B. He), [email protected] (H. Zhou). https://doi.org/10.1016/j.jmgm.2019.107508 1093-3263/© 2019 Published by Elsevier Inc.

survey about their thermodynamic stabilities [15e19]. Among them, it is found that the cagedshaped molecules are more stable than the other forms, such as N18 [16], N20[17,19,] N24 [18,19], N30 [18,19], N3 [19]2, N36 [18,19], N40 [19], N42 [19]. Why big pure nitrogen molecule in cage form is more stable than it in other form, such as cyclic and acyclic forms? The conclusions of Strout’s works [9,18] are that molecules with the most fivedmembered rings (or pentagons) and threedmembered rings (or triangles) in the nitrogen network tend to be the most stable in the cage structure. Strout’s group in their work [18] further came to a conclusion geometrically that the cylinder-shaped (prismatic) nitrogen cages are more stable than that of any other shapes. This conclusion had already been verified by our previous work [19], where we did some comparative studies on the relative stabilities of deferent nitrogen cages (triangles, quadrangles, pentagons and hexagons on both ends) which general formula are (N6)n with an D3d or D3h symmetry, (N8)n with an D4d or D4h symmetry, (N10)n with an D5d or D5h symmetry and (N12)n with an D6d or D6h symmetry, respectively. It is found that the greater the membered rings on the caps, the shorter the prism and the less number the “layer”, and so the less stability of the cages. The optimizations of cages with quadrangles,

2

Q. Guo et al. / Journal of Molecular Graphics and Modelling 96 (2020) 107508

pentagons and hexagons on both ends do not represent corresponding local minima on the potential energy surface at the DFT/ ccdpvDZ computational level when the nitrogen atoms are just more than 36 [19]. In the contrary, the stable configurations of nitrogen cages with triangles (i.e. molecule (N6)n) on both ends can be obtained even if the number of nitrogen atoms reaches 60 [19] or more [20e23]. It should be noted that the stability factors of pure nitrogen molecules aforementioned are speculations based only on the appearances of the geometric configurations. This assumption has some limitations to explain the molecular stability. In our previous work [19,20], intramolecular interactions are verified to be the substantial factor to make the prismatic (N6)n nitrogen cages more stable than the other prismatic nitrogen cages, such as (N8)n, (N10)n and (N12)n, as more dispersion forces (van der Waals forces) were discovered in (N6)n molecules with an D3d or D3h symmetry [19e23] in the same number of nitrogen atoms. So, it can be said that the stability of pure nitrogen cage should depend on which is dominant among intramolecular interactions, ring tensions and dispersion forces (van der Waals forces). Smaller ring tension and larger dispersion force make the nitrogen cage molecules tend to be stable. The tension of the ring is determined by the tension of the chemical bonds that make up the ring,or rather, by the tension of the three single bonds on each nitrogen atom. The three single bonds of nitrogen atom can be regarded as the common edges of the three polygons around it. That the three polygons are formed by the mixture of pentagons and hexagons can minimize the tension of the three bonds [19]. However, in Strout’s work [18], the results show that the orders of stability are N24(D3d) > N24(D4h) > N24(D6d); N30(D3h) > N30 D5h); N36(D6h) simply doesn’t converge and N36(D3d) is a stable structure at the specified computational level, though the order of magnitude of pentagons is opposite to the order of stability above. This is because the dispersion force is found in multiple regions inner N24(D3d), while there is no dispersion force found inner N24(D4h), N24(D6d), N30(D3h), N30(D5h) and N36 (D6h). These findings suggest that in larger nitrogen cage molecules, dispersion forces are the dominant stability factor. This conclusion was verified by our previous work [19e24] in which it can be found that for nitrogen cages with quadrangles, pentagons and hexagons on both ends, the largest stable molecule are only N56(D4h), N30(D5h) and N24(D6d) respectively, and for the nitrogen cages triangles on both ends, the largest stable molecule reported is N90(D3h) or even larger. It is noted that the maximum value of n in the (N6)n molecular sequence with an D3d or D3h symmetry has not yet been determined. In this work, the nitrogen molecules (N6)n (n ¼ 16e35, with D3h or D3d symmetry alternatively) are subjected to computational survey to determine the maximum value of n of (N6)n cages and to probe the mechanism of dispersion forces (van der Waals forces) existing in the stable (N6)n cages. 2. Methods and results For the convenience of the readers, it need to re‒describe the design procedures of molecules (N6)n in our previous work [19], where N6 triangular prism (featured D3h symmetry) was chosen as the parent molecule and a chair‒form six‒membered ring (a hexagonal molecular segments) was chosen as a repeating unit. The parent molecule N6 has two parallel ends, called top end and bottom end respectively. Firstly, the three bonds between the two parallel ends are disconnected, and then one chair hexagon unit is inserted into the parent molecule midway between the two parallel ends. The atoms of the inserted hexagon unit are linked to the atoms in the top end and in the bottom end in an alternating manner. Thus, a new N12 cage with D3d symmetry is formed. In the hexagon unit of N12(D3h) cage, the three atoms linking to the top end atoms are in a same plane, and the remaining three atoms in

the hexagon are another plane, both of which are perpendicular to the C3 principal axis of N12. Analogously, the cage molecule N18(D3h) can be regarded as two chair hexagons unit being inserted symmetrically into the parent N6 molecule. By that analogy, the subsequent molecules N96(D3d), N102(D3h), N108(D3d), N114(D3h), N120(D3d), N126(D3h), N132(D3d), N138(D3h), N144(D3d), N150(D3h), N156(D3d), N162(D3h), N168(D3d), N174(D3h), N180(D3d), N186(D3h), N192(D3d), N198(D3h), N204(D3d) and N210(D3h) can be constructed by inserting 16 to 35 chair hexagons unit into the parent N6 molecule, respectively. For short, the inserted chair hexagons unit is expressed as “Chair” below. A prismatic molecule (N6)n contains n ‒ 1 Chair. From the point of view of molecular construction, every single insertion of a hexagon into parent N6 yields two “atomic planes” which are parallel with each other and perpendicular to the C3 principal axis. For the convenience, the concept of “atomic layer”, or “layer” for short, is introduced to describe these atomic planes. For the original definition of the term “layer”, see reference 19. Since two new layers are formed when one Chair is inserted into the parent molecule, total number of layers is even (a total of 2n) for a prismatic molecule in the (N6)n sequence. Take molecule (N6)16 (D3d) in the (N6)n molecular sequence as example, the Chair (a total of n ‒ 1) and layers (a total of 2  16 ¼ 32) in the molecule and the atoms in each layer are shown in Fig. 1. It is noted that the moiety between layer 1 (top end or left end) and Chair 1, and between Chair n ‒ 1 and layer 2n (bottom end or right end) are both segments of truncated triangular pyramid, called “prismatic cap” in following discussions. The rest is the hexagonal prismatic segment which outside surface is called “prismatic surface”. Neither prismatic cap nor prismatic surface is labelled in Fig. 1 or other figures. On the surface of cage molecules (N6)n optimized at the B3LYP/ cc‒pvDZ level, there are three types of polygons of which two are fused together to form a chemical bonds: triangle, pentagon and hexagon. Triangles are located at the two ends. Six pentagons are laid the side face of two prismatic caps. The remain polygons on the prismatic surface are all boat-form hexagons (Boat for short below) which number is 3(n‒2) derived from the law that there are three Boat around the hexagonal prismatic section (segment) between ith and (i þ 1)th Chair. For example, in optimized N96 (D3d) cage, the three Boat between Chair 8 and Chair 9 are shown in Fig. 2. It is noted that the n ‒ 1 hexagons inserted in N6 keep their chair conformation unchanged in molecule (N6)n. For a three-dimensional molecular graphics of (N6)n which is laid flat, the top view and the front view are identical, both view orientations are perpendicular to the C3 principal axis, while the side view (left view or right view) goes along the C3 principal axis. For all the molecules in (N6)n sequence, the top views and the front views are similar to Fig. 1, which differences only exist in the number of layers or in the length of the prisms. According to the molecular side views (left views or the right views), the molecules in (N6)n sequence are divided into two subsets by their symmetry: One subset is characterized as D3h symmetry and odd number n in the (N6)n sequence. Molecule (N6)17, i.e. N102 (D3h) is the example of this subset (Fig. 3). Another subset is featured with D3d symmetry and even number n in the (N6)n sequence. Molecule (N6)18, i.e. N108 (D3d) is the example of this subset (Fig. 4). The three-dimensional molecular graphics of (N6)n sequence in details can be found in Data-In-Brief (Fig. 123). In this work, the primary simulations software is the wellknown Gaussian09 package [25] which has been widely used in theoretical chemistry research for a long time. The hybrid nonlocal B3LYP density functional theory with Dunning’s correlationconsistent cc‒pvDZ basis set has been carefully picked out for optimizing the structures and vibrations of the (N6)n molecule sequence. In addition, some auxiliary analysis tools, such as natural

Q. Guo et al. / Journal of Molecular Graphics and Modelling 96 (2020) 107508

3

Fig. 1. The front view of the example molecule (N6)16(D3d) optimized at the B3LYP/cc‒pvDZ level.

bond orbital (NBO) analysis [26] and atoms in molecule (AIM) analysis [27], were introduced to obtain the topological properties of bonds in the molecules. The thermodynamic stability is closely related to the topological properties of bonds in a molecule. In AIM analysis, the topological properties of a bond are considered as a topological set of the electron density (r), the Laplacian of r (V2r), the bond critical points (BCPs), the ring critical points (RCPs) and the cage critical points(CCPs), the bond paths (BPs), the ring paths(RPs), the cage paths (CPs) and the distances between the BCPs and the atoms. Wherein the V2r is an important measuring indicator for the charge of the region between two atoms. The V2r > 0 means that the electron density is locally depleted and is typically associated with

Fig. 2. The three boat hexagons around prism segment between Chair 8 and Chair 9.

an interaction (ionic bond, hydrogen bond and van der Waals force) between the internuclear region, whereas the V2r < 0 signifies that the electron density is locally concentrated and characterizes a covalent bond formed in the region between the two nucleus. Also, a parameter “ellipticity” (ε) describing the characteristics of chemical bonds can be derived from the electron density (r), where the greater the ε value is, the more p component the bond has; the less the ε value is, the more s component the bond has. Naturally, a bond intensity can be measured by the r value of its BCP, where the greater the r value is, the more intensive the bond is. For the example molecule (N6)16, the front view of all the critical points and paths as well as the 96 atoms is as Fig. 5. In order to clearly observe the distribution of the two dominating indicators BCPs and BPs in the molecular AIM diagram, the other indicators RCPs, CCPs, RPs and CPs are hidden from view. The BCPs and BPs distribution as well as nucleus of the example molecule (N6)16 viewed from front orientations are shown in Fig. 6. It is noted that the picture elements for large-sized prismatic cages, for example N180 (D3d) (see Fig. 18 in Data-In-Brief.), are too small in the AIM view window as a whole molecule to be labelled and distinguished. This is especially true for even larger molecules, such as N204 (D3d) and N210 (D3h) etc. Originally, the B3LYP functional does not address intermolecular van der Waals interactions well due to its missing long-range correction for the generalized gradient approximation exchange functionals [28]. However, AIM analysis, as a model being used to obtain the electronic density and, hence the dispersion interaction results, is relatively independent of the choice of the functionals as well as the basis sets in the respect of describing the intramolecular interactions [29]. Based on the charge density and the concept of a gradient path, AIM partitions a molecular system into atoms in a natural and rigorous way. Integrating properties over the atomic basins can yield results of the inter- or intra-molecular interactions [30]. This is the reason why the DFT method and the AIM theory were combined in this work to study the weak interactions. 3. Discussion

Fig. 3. Left view of the molecule (N6)

Fig. 4. Left view of the molecule (N6)

17

18

i.e. N102 with the symmetry D3h.

i.e. N108 with the symmetry D3d.

3.1. Geometries optimization and the maximum number n in (N6)n All the geometry structures of prismatic molecules (N6)n (n ¼ 16e35) were optimized at the B3LYP/cc‒pvDZ level (see Fig. 1 of example molecule N96) with the four convergence indexes (Maximum Force, RMS Force, Maximum Displacement, RMS Displacement) are all “YES”. Though further optimizations after vibration calculations show that the “Maximum Displacement” among the four convergence indexes is “NO”, the optimized structures of prismatic molecules (N6)n are taken as reasonable. The vibration mode of the first vibration frequency of all molecules (N6)n is showed as Fig. 7, which indicate that the atoms at both ends of prismatic molecules oscillate in the positive direction of the Yaxis while the atoms in the middle moiety oscillate in the negative direction of the Y-axis. The vibration modes of the first and second vibration frequency is identical. In other words, the macroscopic effect is that the molecules oscillate between linear and curved

4

Q. Guo et al. / Journal of Molecular Graphics and Modelling 96 (2020) 107508

Fig. 5. The 3d side view of the example molecule N96 in AIM window. All the diagram elements of nitrogen atoms (labelled), BCPs (small balls), RCPs (small balls), CCPs (small balls), BPs, RPs, CPs are shown.

Fig. 6. The 3d side view of the example molecule N96 in AIM window. Three diagram elements of nitrogen atoms (labelled balls), BCPs (small balls), and BPs (lines) are shown.

Fig. 7. Schematic diagram of vibration mode of the first vibration frequency of the molecule (N6)

lines. Obviously, linear configuration is stable for all the molecules in (N6)n sequence. The first three lowest vibrational frequencies of the molecules are listed in Table 1. The lowest vibrational frequencies are real at the B3LYP/cc‒pvDZ level except the cage N210, i.e. (N6)35, of which the lowest vibrational frequency value is about 1.44 cm1. The lowest frequencies decrease slowly with the increasing n in (N6)n. A small value of imaginary frequency occurs in (N6)35, which may not necessarily reflect the real vibration situation in the direction of the spatial freedom, because the force of the first two vibration modes of this molecule is very small (about 0.0001 mDyne Å1), and the harmonic frequencies are easily polluted by translations and rotations. Considering the similarity of the first two vibration modes, the first vibration frequency is a small imaginary frequency value (1.44 cm1) in (N6)35, which does not affect the judgment that the molecule is relatively stable at this point. Suppose that the real frequency and the imaginary frequency are token as the criteria for judging whether or not a structure represents a local minimum on potential energy surface, from n ¼ 35 on, the structures of subsequent species (N6)n (n > 35) will not be stable at the B3LYP/cc‒pvDZ level. In fact, many imaginary frequencies occur at the B3LYP/cc‒pvDZ level when n ¼ 36. The lowest vibrational frequency value being about 1.44 cm1 indicates a saddle point in (N6)35 or a metastable structure of (N6)35.

16

i.e. N96 with the symmetry D3d optimized at the B3LYP/cc‒pvDZ level.

Table 1 Total bonds and the first three lowest frequencies of each molecule in (N6)n optimized at the B3LYP/cc‒pvDZ level. n

(N6)n

Freq_1(cm1)

Freq_2(cm1)

Freq_3(cm1)

Bond(n)

16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35

N96 N102 N108 N114 N120 N126 N132 N138 N144 N150 N156 N162 N168 N174 N180 N186 N192 N198 N204 N210

19.3173 18.3806 16.4253 14.9928 13.7555 12.4984 11.3334 10.2464 9.6039 8.5088 8.4170 7.3657 6.8679 6.7469 5.8848 4.0850 4.4669 3.9579 2.7756 1.4414

19.8090 18.8016 16.8996 15.4424 14.1378 12.7753 11.6149 10.6119 9.9744 8.9552 8.8715 7.9548 7.6793 7.3883 6.9060 6.5891 5.9551 5.2682 5.0142 4.7971

27.2230 33.6139 30.5817 34.8133 32.0363 31.4594 28.7725 26.3435 24.4696 16.8424 21.1947 19.4250 18.1584 17.1544 15.8292 14.0161 13.0338 12.5047 11.3819 10.2565

144 153 162 171 180 189 198 207 216 225 234 243 252 261 270 279 288 297 306 315

Q. Guo et al. / Journal of Molecular Graphics and Modelling 96 (2020) 107508

So, n ¼ 35 can be regarded as the threshold value of stable structure in (N6)n molecular sequence. Furthermore, all the 20 optimized species (N6)n maintain their prismatic structures (see Fig. 1 above and Fig. 123 in Data-In-Brief). Each nitrogen atom has three single bonds linking to its three adjacent N atoms. The number (column M(n) in Table 1) of covalent bond in (N6)n molecule follows formula (F1):

BondðnÞ ¼


 3 6  n ¼9n 2

(F1)

The topological structure and number Bond(n) of covalent bonds in (N6)n molecule can further be confirmed by NBO and AIM analyses below. Based on the discussions above, a conclusion that the maximum number n in (N6)n sequence, of which structure is stable at the B3LYP/cc‒pvDZ level, is 35 can be drawn.

5

reasonable, including the largest sized molecule (N6)35, i.e. N210. For the convenience of categorizing the vast amount of data yielded from the structure optimizations, bonding types were defined in the same way of our previous works [19e24]: Each bond can be considered as a shared edge of two adjacent polygons on the prismatic caps and prismatic surface of a molecule. The bonding properties are related to the shared edge of the two polygons. If the edge numbers of two adjacent polygons are m (m ¼ 3e6) and n (n ¼ 3e6), respectively, the shared edge of the two polygons is defined as the mn (m ≤ n) bond, labelled mnB, in this work. In each molecule of (N6)n (n ¼ 16e35) sequence, there are five different types of bonds which belong to five sets, called mnB set. As an example, the five sets of bonds in molecule (N6)16, i.e. N96, are categorized in Table 3. 3.3. NBO and AIM analyses

3.2. Energies The various thermal energies, such as energies, enthalpies and free energies, of the stable molecules, including the parent molecule N6 are listed in Table 2. DEi represents the differentials between the molecular thermal energy and its previous one. The energy of a molecule in (N6)n sequence decreases linearly by a certain difference, DEi ¼ 328.155674 to 328.156608 au. According to the rule of molecular construction aforementioned, one Chair is added to each successive molecule in the (N6)n molecular sequence. So, DEi is the energy contribution of one Chair fragment to the successive molecule. EN6 denotes the thermal energy of parent molecule N6. The highest value of DEi (328.155674 au) is still lower than that of EN6 (328.092092 au). DDEi is defined to indicate this differential. That is, suppose a thermochemical reaction: (N6)i-1 þ N6 (Chair) / (N6)i, and keep the shapes of (N6)i-1 and N6 (Chair) unchanged in (N6)i

DDEi amount of heats would be released since all the values of them are negative. DDEi also can be taken as the stabilization energy from (N6)i-1 to (N6)i. The existence of stabilization energy makes the stability of molecule in (N6)n molecular sequence

NBO analysis can yield information about molecular orbitals which include the bond and anti-bond orbitals, and lone‒pair and Rydberg orbitals, etc. NBO results show that there are 6n lone‒pair orbitals distributed over each atom in (N6)n cage, i.e. each atom has only one lone‒pair orbital filled with two electrons; there are three bond orbitals on each atom, which orbitals overlap with three adjacent N atoms to form three sN‒N single bonds. That is, the bond number in each (N6)n cage from NBO analysis also follows the same formula (F1) as it in geometric optimization (see Table 4). In other words, the bonding orbitals obtained from the NBO analysis are matched with the chemical bonds obtained from the geometric structure optimization. Total number of banding orbitals for the twenty (N6)n (n ¼ 16e35) cages can be calculated by formula (F2) Bond

Total

¼

35 X

MðnÞ ¼

n¼16

35 X

! 9n

¼ 4590

(F2)

n¼16

BCPs from AIM analysis fall into two broad categories. One is the BCPs with V2r < 0 which indicate covalent bonds. The other is the BCPs with V2r > 0 which imply ionic bond, hydrogen bond, dispersion and van der Waals forces. After classifying and sorting, it is found that the number of BCPs with V2r < 0 in each (N6)n cage is also calculated by formula (F1). Total number of BCPs with V2r < 0

Table 2 Various energies of each molecule in (N6)n and the parent molecule N6 optimized at the B3LYP/cc‒pvDZ level. n

(N6)n

Energy(au)

Enthalpy(au)

Free Energy(au)

DEi (au)

DDEi (au)

DDEi (KJ mol-1)

1 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35

N6 N96 N102 N108 N114 N120 N126 N132 N138 N144 N150 N156 N162 N168 N174 N180 N186 N192 N198 N204 N210

328.092092 5250.577390 5578.733496 5906.889443 6235.045462 6563.201390 6891.357335 7219.513240 7547.669204 7875.825150 8203.980846 8532.136943 8860.292780 9188.448760 9516.604951 9844.760738 10172.916974 10501.072648 10829.228934 11157.384779 11485.541387

328.091148 5250.576446 5578.732552 5906.888499 6235.044517 6563.200446 6891.356390 7219.512295 7547.668259 7875.824206 8203.979902 8532.135999 8860.291836 9188.447816 9516.604007 9844.759793 10172.916030 10501.071704 10829.227990 11157.383835 11485.540442

328.120904 5250.703362 5578.864982 5907.027369 6235.189268 6563.351587 6891.513799 7219.676212 7547.838432 7876.000683 8204.163901 8532.325494 8860.488073 9188.650374 9516.812017 9844.974888 10173.137387 10501.300008 10829.462039 11157.624652 11485.783234

328.156106 328.155947 328.156019 328.155928 328.155945 328.155905 328.155964 328.155946 328.155696 328.156097 328.155837 328.155980 328.156191 328.155787 328.156236 328.155674 328.156286 328.155845 328.156608

0.064014 0.063855 0.063927 0.063836 0.063853 0.063813 0.063872 0.063854 0.063604 0.064005 0.063745 0.063888 0.064099 0.063695 0.064144 0.063582 0.064194 0.063753 0.064516

168.07 167.65 167.84 167.60 167.65 167.54 167.70 167.65 166.99 168.05 167.36 167.74 168.29 167.23 168.41 166.93 168.54 167.38 169.39

DEi ¼ Ei e Ei-1. DDEi ¼ DEi e EN6.

6

Q. Guo et al. / Journal of Molecular Graphics and Modelling 96 (2020) 107508 Table 3 Subtotal results of covalent bonds, lone pairs and dispersion interactions in the example N96 cage (D3d symmetry) obtained from geometry optimization, NBO and AIM analyses at the B3LYP/cc‒pvDZ level, respectively. Type

Number

Bond

35B 55B 56B

6 6 12

66B1

42

66B2

78

V1* V2

6 36

LP1 LP2 LP3 LP4

6 6 6 78

N1eN2, N1eN3, N2eN3, N94eN95, N94eN96, N95eN96 N1eN4, N2eN6, N3eN8, N89eN94, N91eN95, N93eN96 N4eN5, N4eN9, N5eN6, N6eN7, N7eN8, N8eN9, N88eN89, N88eN93, N89eN90, N90eN91, N91eN92, N92eN93 N5eN11, N7eN13, N9eN15, N10eN16, N12eN18, N14eN20, N17eN23, N19eN25, N21eN27, N22eN28, N24eN30, N26eN32, N29eN35, N31eN37, N33eN39, N34eN40, N36eN42, N38eN44, N41eN47, N43eN49, N45eN51, N46eN52, N48eN54, N50eN56, N53eN58, N55eN60, N57eN62, N59eN64, N61eN66, N63eN68, N65eN70, N67eN72, N69eN74, N71eN76, N73eN78, N75eN80, N77eN82, N79eN84, N81eN86 N83eN88, N85eN90, N87eN92 N10eN11, N10eN15, N11eN12, N12eN13, N13eN14, N14eN15, N16eN17, N16eN21, N17eN18, N18eN19, N19eN20, N20eN21, N22eN23, N22eN27, N23eN24, N24eN25, N25eN26, N26eN27, N28eN29, N28eN33, N29eN30, N30eN31, N31eN32, N32eN33, N34eN35, N34eN39, N35eN36, N36eN37, N37eN38, N38eN39, N40eN41, N40eN45, N41eN42, N42eN43, N43eN44, N44eN45, N46eN47, N46eN51, N47eN48, N48eN49, N49eN50, N50eN51, N52eN53, N52eN57, N53eN54, N54eN55, N55eN56, N56eN57, N58eN59, N58eN63, N59eN60, N60eN61, N61eN62, N62eN63, N64eN65, N64eN69, N65eN66, N66eN67, N67eN68, N68eN69, N70eN71, N70eN75, N71eN72, N72eN73, N73eN74, N74eN75, N76eN77, N76eN81, N77eN78, N78eN79, N79eN80, N80eN81, N82eN83, N82eN87, N83eN84, N84eN85, N85eN86, N86eN87 N4╍N10, N6╍N12, N8╍N14, N82╍N93, N84╍N89, N86╍N91 N11╍N17, N13╍N19, N15╍N21, N16╍N22, N18╍N24, N20╍N26, N23╍N29, N25╍N31, N27╍N33, N28╍N34, N30╍N36, N32╍N38, N35╍N41, N37eN43, N39╍N45, N40╍N46, N42╍N48, N44╍N50, N47╍N53, N49╍N55, N51╍N57, N52╍N63, N54╍N59, N56╍N61, N58╍N69, N60╍N65, N62╍N67, N64╍N75, N66╍N71, N68╍N73, N70╍N81, N72╍N77, N74╍N79, N76╍N87, N78╍N83, N80╍N85 N1, N2, N3, N94, N95, N96 N4, N6, N8, N89, N91, N93 N5, N7, N9, N88, N90, N92 N10 ~ N87

* The intensity dispersion interactions V1 is weaker than that of V2. Note: mnB (m ¼ 3, 5, 6, and n ¼ 5, 6) denotes covalent bonds which come from geometry optimization, NBO and AIM analyses at the B3LYP/ cc‒pvDZ level; 66B1 and 66B2 are two types of 66B bonds; V1 and V2 denote van der Waals force or dispersion interactions which come only from AIM analysis; LPn (n ¼ 1e4) denotes lone pairs which come only from NBO analysis.

Table 4 The number of covalent bonds and their classification in NBO and AIM Analyses, and the number of van der Waals forces in AIM analysis at the B3LYP/cc‒pvDZ level. n

(N6)n

Bond(n)

35B

55B

56B

66B1

66B2

vdW(n)

16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35

N96 N102 N108 N114 N120 N126 N132 N138 N144 N150 N156 N162 N168 N174 N180 N186 N192 N198 N204 N210

144 153 162 171 180 189 198 207 216 225 234 243 252 261 270 279 288 297 306 315

6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6

6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6

12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12

42 45 48 51 54 57 60 63 66 69 72 75 78 81 84 87 90 93 96 99

78 84 90 96 102 108 114 120 126 132 138 144 150 156 162 168 174 180 186 192

42 45 48 51 54 57 60 63 66 69 72 75 78 81 84 87 90 93 96 99

for the twenty (N6)n (n ¼ 16e35) cages can be gotten by formula (F2), i.e. 4590. In a word, the connections of atoms (bonds) obtained from geometric optimization, bonding orbitals of NBO analysis and BCPs with V2r < 0 of AIM analysis, respectively, are

entirely consistent. So, these three terms can be combined into one, collectively referred to as bond simply in this paper. The number of BCPs with V2r > 0 in each of (N6)n cage can be obtained by formula (F3)

vdWðnÞ ¼ 3  ðn  2Þ

(F3)

and the total number of BCPs with V2r > 0 for the twenty (N6)n (n ¼ 16e35) cages can be calculated by formula (F4)

TotalvdW ¼

35 X n¼16

vdWðnÞ ¼

35 X

3 ðn  2Þ ¼ 1410

(F4)

n¼16

The number of chemical bonds which are determined by geometric optimization, NBO and AIM analyses, and their classification, as well as the number of van der Waals forces in each molecule are listed in Table 4. The most important topological parameters related to chemical bonds are those of bond lengths obtained from structural optimization, bonding orbital occupancies and bonding orbital energies obtained from NBO analysis, BP lengths, ellipticity (ε) values, electron density (r) and Laplace value (V2r) obtained from AIM analysis. If the data of bonds (from optimization), bond orbitals (from NBO) and BCPs (from AIM) were listed separately, there would be 4,590  3 ¼ 13,770 data lines besides dozens of the header lines. For each molecule, the topological parameters of the bonds in

Q. Guo et al. / Journal of Molecular Graphics and Modelling 96 (2020) 107508

geometric optimization, bonding orbitals in NBO analysis and BCPs in AIM analysis can be integrated into one table due to their consistency. Even so, in addition to the 1410 van der Waals forces in the twenty molecules, there are still as many as 4,590 þ 1410 ¼ 6000 data rows in twenty tables(see Data-In-Brief). In the text, the properties and functions of chemical bonds are discussed based on the example molecule N96 and representative bonds of bond categories: 35B, 55B, 56B, 66B1, 66B2 and vdW since the bonding properties of the same bond sets either in the same molecule or in the twenty different molecules are almost identical. It is impossible for the one-to-one comparisons of the bonding parameters among the sequence molecules because the number of bonds varies from molecule to molecule. However, the number of bond sets in the sequence molecule is the same. In order to verify the similarity possessed by bond topological parameters of the same bond set among the sequence molecules, the topological parameters of example molecule N96, which is also the smallest sized molecule in the sequence, and the corresponding bond topological parameters (in bold and italic) of the largest sized molecule N210, are comparably listed in Table 5. It is found that the topological parameters of the same bond sets in N96 and N210 have almost identical values. Thus, a conclusion that other molecules in the sequence have the same type of bond topological parameters as N96, can be logically deduced. It is therefore perfectly reasonable to use the representative bonds of the example molecule N96 to discuss the bonding properties of the whole (N6)n (n ¼ 16e35) sequence. There are six 35B bonds which are located at the two end layers of the prism in each molecule of the (N6)n sequence,. N1eN2 is chosen as the representative bond of this set. Similarly, six 55B bonds lie between the first two layers (layers 1 and 2) and between the last two layers (layers 2n‒1 and 2n) (see Fig. 1), respectively. N1eN4 is chosen as the representative bond. Twelve 56B bonds lie in the first chair hexagon unit inserted (Chair 1) and the last chair hexagon unit inserted (Chair n‒1), respectively. N4eN5 is a sample. The six 66B1 (denoted as 66B1a) bonds between layers 3 and layer 4 and between layer 2n‒3 and layer 2n‒2, taking N5eN11 as example, differ in nature from the other 66B1 (denoted as 66B1b) bonds, taking N46eN52 as sample. Although these six 66B1a bonds are geometrically classified as 66B1 bonds, they’re more like 56B bonds. All 3(n‒2) 66B1 bonds are parallel to the C3 principal axis (see Fig. 1). The largest number (6(n‒3)) of bonds in a molecule is 66B2 which are the sides of chair hexagons inserted except those of Chair 1 and Chair n‒1. The most typical 66B2 bonds is N46eN47 which is chosen as the sample of this set. Similar to the 66B1 bonds,

7

total 3(n‒2) van der Waals forces (vdW) in the molecule are divided into two subsets: one is v1 subset, represented by N4╍N10, of which six van der Waals interactions are located between Chair 1 and Chair 2, and between Chair n‒ ‒2 and Chair n‒1, respectively; the remain 3(n‒2)‒6 van der Waals interactions are classified as v2 subset, represented by N40╍N46. It is noted that an BP length is the sum of the distance between an BCP and its two adjacent atoms instead of the distance between two adjacent atoms. BCP tends to deviate from connection line between two atoms (ε s 0, i.e. BP has some curvature). For an BCP in the prismatic molecule, the BP length is about 0.002 Å larger than the corresponding bond length. Such a small difference is negligible. Therefore, an BP length can be regarded as the approximate bond length, as each of BP looks like quasi beeline. The distance between two atoms refers to the length of the covalent bond for V2r < 0.0, and of the BP for V2r > 0.0 in this work. About the six 35B bonds in the example molecule (N6)16, the bond length of sample bond N1eN2 is the longest (1.529 Å, see Table 5), its orbital energy (OE) is the highest (0.70678 au) and the electron occupancy of the orbital is the lowest (1.954), which means 2.0e1.954 ¼ 0.046 electrons are delocalized in the molecule. The low electron occupancy means that the electrons of this bond as a donor are transferred to the acceptors. For donor/acceptor pair with stabilization energy E(2) > 1.0 kcal mol1, 1.43 kcal mol1 for sN1‒N2/s*N5eN11 and 1.34 kcal mol1 for both sN1‒N2/s*N1eN3 and sN1‒N2/s*N2eN3, indicate that a certain amount of electrons on sN1‒N2 orbital trend to be transferred to s*N5eN11, s*N1eN3 and s*N2eN3 orbitals. Or say the N1eN2 bond has an interaction with N5eN11, N1eN3 and N2eN3 bonds. These interactions contribute to the stability of the molecule. Besides, the electron transfer of lone pairs on both atoms can also contribute to stability of the molecule. Lone pair favors to transfer electron to anti-bonds formed by its three ortho position atoms, take the lone pair on N1 as example, the ortho‒position atoms are N2, N3 and N4, and the anti-bond orbitals are s*N2eN6, s*N3eN8, s*N4eN5 and s*N4eN9. The E(2) of LP / s*NeN are about 2.1 kcal mol1 (see rows of orbital number 241 in Table 6). For donor/acceptor pair with stabilization energy E(2) < 1.0 kcal mol1, the interactions can be ignored. Though the E(2) of sN1‒N2/RY* N3 (anti-Rydberg single center bond) donor/ acceptor pair is greater than 1.0 kcal mol1, the interaction between them can be ignored because the existence of a big energy gap (E(j) e E(i) ¼ 2.02 au) lead to electron transferring with great difficulty. So, the other s / Rydberg* donor/acceptor pairs are no longer listed in Table 6. The s attribute of this covalent bond is confirmed in the NBO and AIM analyses. However, the bond bending in NBO

Table 5 Covalent bond lengths, orbital occupancy energies OE(au), lengths of bond path, ellipticities (ε), electron densities (r) and Laplace of r (V2r) in the N96 and N210 (in BOLD and ITALIC, as a comparison with those of N96) cages at the B3LYP/cc‒pvDZ level.

35B (N1eN2) 55B (N1eN4) 56B (N4eN5) 66B1a (N5eN11) 66B1b (N46eN52) 66B2 (N46eN47) V1(N4╍N10) V2(40╍46)

N96(D3d) N210(D3h) N96(D3d) N210(D3h) N96(D3d) N210(D3h) N96(D3d) N210(D3h) N96(D3d) N210(D3h) N96(D3d) N210(D3h) N96(D3d) N210(D3h) N96(D3d) N210(D3h)

Length (Å)

Occup.

OE(au)

BP Length(Å)

ε

r(e/Bohr3)

V2r

1.529 1.529 1.428 1.427 1.483 1.483 1.483 1.483 1.491 1.491 1.463 1.463

1.954 1.954 1.974 1.974 1.976 1.976 1.970 1.970 1.966 1.966 1.984 1.984

0.70678 0.70665 0.89777 0.89832 0.84599 0.84625 0.86759 0.86759 0.86432 0.86467 0.88118 0.88193

1.532 1.532 1.429 1.429 1.485 1.484 1.485 1.485 1.492 1.492 1.464 1.464 2.418 2.418 2.336 2.336

0.182 0.183 0.105 0.105 0.106 0.107 0.199 0.199 0.208 0.208 0.091 0.091 3.233 3.223 0.049 0.049

0.247 0.247 0.323 0.323 0.285 0.285 0.282 0.282 0.278 0.278 0.298 0.298 0.024 0.024 0.03 0.03

0.247 0.247 0.617 0.618 0.456 0.457 0.498 0.498 0.494 0.494 0.515 0.515 0.143 0.143 0.162 0.161

8

Q. Guo et al. / Journal of Molecular Graphics and Modelling 96 (2020) 107508

Table 6 The second order perturbation theory analysis of Fock matrix in NBO for selected orbitals of example molecule N96. E(2) is the stability energy and E(j)-E(i) is the energy gap between acceptor and donor. Orb. Num.

Donor (i)

Acceptor*(j)

E(2) (kcal mol1)

E(j)-E(i) (au)

1

N1eN2

RY* N3 N1eN3, N2eN3 N1eN4, N2eN6 N4eN5, N5eN6 N5eN11 N5eN11, N9eN15 N1eN2 N11eN12 N1eN4, N2eN6 N10eN16, N12eN18 N29eN35, N33eN39 N41eN47, N45eN51 N34eN39, N47eN48 N34eN35, N50eN51 N41eN42, N52eN57 N41eN47, N45eN51 N53eN58, N57eN62 N2eN6, N3eN8 N4eN5, N4eN9 N1eN2, N1eN3 N5eN6, N8eN9 N4eN9, N6eN7 N10eN11, N11eN12 N34eN35, N34eN39 N41eN42, N44eN45 N46eN52 N34eN40 N47eN48, N50eN51 N52eN53, N52eN57

1.05 1.34 0.73 0.87 1.43 2.65 1.11 0.63 1.05 2.23 2.10 2.11 0.59 0.59 0.59 2.10 2.11 2.11 2.04 3.87 5.32 7.80 2.36 2.27 8.17 1.63 1.63 8.16 2.28

2.02 0.74 0.87 0.80 0.81 1.00 0.88 0.96 1.03 0.96 0.96 0.96 0.99 0.99 0.99 0.96 0.96 0.74 0.67 0.53 0.59 0.58 0.61 0.59 0.59 0.58 0.58 0.59 0.59

3 7

N1eN4 N4eN5

10

N5eN11

54

N34eN40

61 62 70 72

N40eN41 N40eN45 N46eN47 N46eN52

241

N1

244

N4

245

N5

280

N40

286

N46

analysis shows a large angle deviation (19.6 for both hybrid orbitals of the two atoms N1 and N2, see row of orbital number 1 or 144 in Table 7) from the hybrid center (marked with q and f), while the angle deviation is very close to 0 for the localized s single bond, for example, no angle deviation value is given for hybrid orbital of N5 in bonding orbital N5eN11 (see row of orbital number 7 in Table 7) because the printing threshold value of angle deviation is 0.5 . The negative V2r value (0.247) and the larger ε value (0.182) of its BCP’s should be regarded as partial p-orbital components are superimposed on the sN1N2 single bond to strengthen it. This view is supported by a larger r value (0.248 e/Bohr3) of the region. The longest bond length and the highest OE of this bond may be due to the larger tensions of bonds. All the nitrogen atoms in the molecule are sp3 hybridized, which means that the most appropriate bond angle should be close to the angle of 109.5 in tetrahedral configuration. For atom N1, the bond angles are :N2eN1eN3 z 60 ,:N2eN1eN4 ¼ :N3eN1eN4 z 106.4 , respectively. Obviously, the angle :N2eN1eN3z 60 , i.e. the angle between N1eN2 and N1eN3 both being 35B, is far from appropriate bond angle of 109.5 . This extremely small angle causes a lot of small-angle strain (repulsion) between the 35Bs. To alleviate this tension, lengthening bond length (1.484 Å for typical sN1N2 single bond) is one of the effective ways. The presence of tension can weaken the bond and consequently lower the relative stability of the molecule. On the other hand, the double bond tendency enhances the bond intensity and increases the stability of the molecule. It can be speculated that, between these two counter stabilizing factors, the double bond tendency is the dominant one. Also, the N1eN2 bond is a shared edge between a triangle and a pentagon. The triangle acting as stabilizing factor [18] in nitrogen cages may be caused by the double bond tendency. N1eN4 is the sample of the six 55Bs in the example molecule N96. The length is the shortest (1.428 Å) and its OE is the lowest

Table 7 Bond bending (deviations from line of nuclear centers) measured by angular (q, f in  ) and deviations (dev in  ) between hybridized orbitals and nuclear centers of the sample bond in the example molecule N96. No.

Bond

1 3 7 10 70 72 144

N1eN2 N1eN4 N4eN5 N5eN11 N46eN47 N46eN52 N95eN96

Line of Center

Hybrid Orbital of Atom 1

Hybrid Orbital of Atom 2

q

f

q

f

dev

q

f

dev

98.1 97.0 127.6 89.2 126.0 90.0 81.9

270.0 17.8 296.6 359.9 249.3 0.0 90.0

117.7 97.9 126.5 – 124.7 90.5 62.3

268.8 20.2 301.7 – 249.4 1.3 88.8

19.6 2.6 4.2 – 1.3 1.5 19.6

101.4 84.2 52.4 92.1 52.8 90.6 78.6

91.1 194.6 112.1 179.7 68.7 178.6 271.1

19.6 3.4 3.5 1.3 1.3 1.5 19.6

(0.89777 au) in the molecule. Its electron occupancy of the orbital is 1.974, which indicates that the number of nonlocalized electrons is small (2.0e1.974 ¼ 0.026). The stabilization energy E(2) is 2.65 kcal mol1 for both sN1‒N4/s*N5eN11 and sN1‒N4/s*N9eN15 (see rows of orbital number 3 in Table 6), indicating that there are interactions between N1eN4 and its two meta‒position bonds N5eN11 and N9eN15. The bond is confirmed as a covalent bond by the NBO analysis and by the AIM analysis as the BCP’s V2r is negative (0.617). The electron density r value (0.323 e/Bohr3) is the greatest among all of the BCPs’ electron densities, which indicates that the BCP has the highest electron density and, consequently, the corresponding covalent bond should be the strongest in the molecule. The bond angles, formed by N1eN4 and its four adjacent bonds, are :N1eN4eN5 ¼ :N1eN4eN9 z 101.2 ,:N2eN1eN4 ¼ :N3eN1eN4 z 106.4 , respectively. Both angle values are not too far from appropriate tetrahedral bond angle of 109.5 , indicating that the tensions (repulsion) between N1eN4 and its four adjacent bonds are small. The small tensions between N1eN4 and its four adjacent bonds can be verified by the bond bending of N1eN4, which angle deviations are only 2.6 for hybrid orbitals of atom N4 and 3.4 for hybrid orbitals of atom N5 (see row of orbital number 3 in Table 7) from the hybrid center. In addition, the bond length of N1eN4 is only 1.428 Å, which is between a normal resonant bond NN (1.316 Å) and a normal sNN single bond (1.484 Å). By the aforementioned definition, the number of nonlocalized electrons of N4 is the greatest (2.0e1.893 ¼ 0.107). In addition, the lone pair on N4 besides N1 has relatively strong interactions with s*N1eN2 (E(2) ¼ 3.87 kcal mol1), s*N1eN3 (E(2) ¼ 3.87 kcal mol1), s*N5eN6 (E(2) ¼ 5.32 kcal mol1), s*N8eN9 (E(2) ¼ 5.32 kcal mol1) (see rows of orbital number 244 in Table 6). The nonlocalized electrons strengthen the bonds between atoms N4 and the neighboring atoms. Furthermore, this bond is a shared edge between two adjacent pentagons. The pentagon being a dominant stabilizing factor [18] in nitrogen cages may be due to the increased nonlocalized electrons of the lone pairs on the atoms in the pentagon. There are twelve 56B bonds (N4eN5 as the sample), which are the sides of Chair 1 and Chair n‒1 in each molecule of (N6)n sequence. Bond length of the sample N4eN5 is 1.483 Å, which length is much closer to the normal sNN single bond (1.484 Å). Its lower OE (0.84599 au) higher electron occupancy (1.976) (see Table 5) suggest that the number of nonlocalized electrons is small (2.0e1.976 ¼ 0.024). The stabilization energy E(2) is 1.11 kcal mol1 for sN4‒N5/s*N1eN2 and only 0.63 kcal mol1 for sN4‒ N5/s*N11eN12 (see rows of orbital number 3 in Table 6), indicating that the interactions between N4eN5 and its two meta‒position bonds N1eN2 and N11eN12 are weak. The bond is confirmed as a covalent bond by the NBO analysis and by the AIM analysis as the BCP’s V2r is negative (0.456). The electron density r value (0.285

Q. Guo et al. / Journal of Molecular Graphics and Modelling 96 (2020) 107508

e/Bohr3) indicates that the corresponding covalent bond should be relatively strong. The three bond-angles :N1eN4eN5 ¼ :N1eN4eN9 z 101.2 , :N5eN4eN9 z 109.3 with vertex N4 and the three bond-angles:N4eN5eN6 z 103.5 , :N4eN5eN11 ¼ :N6eN5eN11 z 112.2 with vertex N5 are close to the appropriate tetrahedral bond angle of 109.5 , implying that the tensions (repulsive forces) between N4eN5 and its four adjacent bonds are small. The small repulsions between N4eN5 and its four adjacent bonds can be verified by low ε value (0.106) and by the bond bending of N4eN5 in NBO analysis, which angle deviations are 4.2 for hybrid orbitals of atom N4 and 3.5 for hybrid orbitals of atom N5 (see row of orbital number 7 in Table 7) from the hybrid center. N4eN5 is a typical sN‒N single bond. Its low orbital energy and low repulsions with adjacent bond contribute to the stability of the whole molecule. 66B1a are very similar to the twelve 56B bonds (N4eN5 as the sample bond) because their bond lengths, electron occupancies, ellipticity values (ε) and electron density values (r) are similar. The differences are those: i) the OE (0.86759 au) of N5eN11 is lower than that of N4eN5 (0.84599 au); ii) The stabilization energy E(2) is 1.05 kcal mol1 for both sN5‒N11/s*N1eN4 and sN5‒ 1 for both sN5‒N11/s*N10eN16 N11/s*N2eN6, and 2.23 kcal mol and sN5‒N11/s*N12eN18 (see rows of orbital number 10 in Table 6); iii) The three bond-angles:N4eN5eN6 z 103.5 , :N4eN5eN11 ¼ :N6eN5eN11 z 112.2 with vertex N5 and the three bond-angles :N5eN11eN10 ¼ :N5eN11eN12 z 110.5 , :N10eN11eN12 z 108.6 with vertex N11 are very close to the appropriate angle of 109.5 ; iv) The bond bending of N5eN11 in NBO analysis show that no angle deviations for hybrid orbitals of atom N5 and only 1.3 for hybrid orbitals of atom N11 (see row of orbital number 10 in Table 7) from the hybrid center. All the points mentioned above imply that the molecular orbital of sN5‒N11 is better overlapped by the sp3 hybrid orbitals of atoms N5 and N11 and contributes to the stability of the molecule. In molecule (N6)n, there are totally 3  (n‒2) 66B1 bonds those are parallel to the C3 principal axis. Except the six 66B1a bonds, the remain 3  (n ‒ 2) ‒ 6 or 3  (n ‒ 4) bonds are 66B1b by getting rid of the six 66B1a bonds from totally 3  (n-2) 66B1 bonds. The differences between 66B1b and 66B1a are the bond lengths (1.491 Å for N46eN52 and 1.483 Å for N5eN11) and the electron occupancies (1.96564 for N46eN52 and 1.96976 for N5eN11), the other topological parameters, such as OE, ε, r, V2r, E(2) and angle deviations are similar greatly. In principal, all the 66B1 (i.e. 66B1a and 66B1b) and 56B bonds have the equivalent contribution to the stability of the molecule. In each molecule of (N6)n sequence, the largest number of bonds is 66B2 (N46eN47 is the sample bond) which is counted by 6  (n‒ 3). 66B2 are the edges of inserted Chair except that the first and the last inserted hexagon unit, i.e. Chair 1 and Chair n‒ ‒1. The bonding orbital of N46eN47 has high electron occupancies (1.984) (Table 5), low ε and small angle deviations (1.3 ) for hybrid orbitals of both atoms N46 and N47 (see row of orbital number 70 in Table 7) from the hybrid center, and low stability energy E(2), implying this bond is highly localized. Its shorter bond length (1.463 Å), low orbital energy OE (0.88118 au) and high electron density r value (0.298 e/ Bohr3) suggest a considerably strong bond and subsequently contribute to the stability of the molecule.

9

forces located between Chair 1 and Chair 2, and between Chair n‒2 and Chair n‒1 are labelled in v1(N4╍N10 as the sample). The remain vdWs are labelled as v2 (N40╍N46 as the sample) in some tables and figures. V2r > 0 (V2r ¼ 0.143 for N4╍N10 and V2r ¼ 0.162 for N40╍N46) indicates that they are both weak interactions (van der Waals interactions or dispersion forces). The BP length (like the bond length of a covalent bond) of N4╍N10 (2.418 Å) being longer than that of N40╍N46 (2.336 Å), and the r value of N4╍N10 (0.024 e/Bohr3) being less than that of N40╍N46 (0.030 e/Bohr3), imply that the interaction of N4╍N10 (located in the ends) is weaker than that of N40╍N46 (located in the middle part). The van der Waals forces or dispersion interactions only appear in the boat hexagons between the two para‒position (i.e. positions 10 and 40 of hexagon) nitrogen atoms which have the proper orientations to effect the interaction (see Figs. 5 and 6). Each boat hexagon has only one vdW. In the molecules, each nitrogen atom is sp3 hybridized and four sp3 hybrid orbitals formed. Among the four sp3 hybrid orbitals, three single electron sp3 hybrid orbitals of a nitrogen atom are paired with the single electron sp3 orbitals of the three surrounding nitrogen atoms to form three covalent bonds. The remain one sp3 is occupied by lone pair electrons. Taking the boat formed by atoms N40, N41, N47, N46, N51, N45 (Fig. 8) as example, three single electron sp3 hybrid orbitals of atom N40 are paired with the single electron sp3 orbitals of the atoms N34, N41, N45, respectively, and form N34eN40, N40eN41, N40eN45. The lone pair orbital of N40 is denoted graphically as P1 in Fig. 8. Just like the atom N40, the lone pair orbital of N46 is denoted graphically as P2 in Fig. 8. The atoms N40 and N46 are located in positions 10 and 40 of the boat hexagon, respectively. The six orbital angles (including lone pair orbital P1) of atom N40 and the six orbital angles (including lone pair orbital P2) of atom N46 are listed in Table 8. According to the angles between orbitals, orbitals N34eN40, P1, P2, N46eN52 are coplanar. From a geometric point of view, orbitals P1 and P2 tend to be close to each other. The proximity effect of nitrogen atom at positions 10 and 4’ of the boat hexagon may cause the intra-molecular interaction (van der Waals force or dispersion force). Suppose P1 is the orbital pivot of lone pair, which is situated at two-thirds of a standard single NeN bond from N40 to P1, the distance between P1 and P2 can be calculated proximately as: BP length of N40╍N46 ¼ 2.336 Å, NeN bond length ¼ 1.484 Å, :N34eN40eP1 ¼ :N52eN46eP2z 110.2 , LP1‒P2 ¼ 2.336e1.484  2  2  |cos(110.2)| ÷ 3 ¼ 1.653 (Å) LP1‒P2 can be regarded as the interaction distance between orbitals P1 and P2. It is larger than the length of standard NeN single bond (1.484 Å), but smaller than the maximum weak interaction

3.4. van der waals interactions in the molecules What’s most exciting is the finding of the same amount of van der Waals forces (vdW, mostly the van der Waals is the intramolecular dispersion interaction) as 66B1 (totally 3  (n‒2)) in AIM analysis of this molecule sequence. Just like the 66B1 bonds, there are two kinds of vdW in the molecule: the van der Waals

Fig. 8. The sketch of the lone pair sp3 orbitals on atoms N40 and N46. A and B are the projection point of the overlap center of the hybrid orbital on line N34eN40╍N46eN52.

10

Q. Guo et al. / Journal of Molecular Graphics and Modelling 96 (2020) 107508

Table 8 sp3 hybrid orbital angles (or bond angles) of atoms N40 and N46 in the example molecule N96. P1 and P2 are the lone pairs on atoms N40 and N46, respectively.

Table 9 The six warp yarns (edges of prismatic surface) formed alternately by covalent bonds () and van der Waals forces (╍) on the side face of the example molecule N96.

Bond angle

Degree ( )

Bond angle

Degree ( )

Number

Warp Yarns (Edges of Prismatic Surface)

:N41eN40eP1 :N45eN40eP1 :N34eN40eP1 :N34eN40eN41 :N34eN40eN45 :N41eN40‒45

110.2 110.2 105.1 111.6 111.6 108.1

:N47eN46eP2 :N51eN46eP2 :N52eN46eP2 :N47eN46eN51 :N47eN46eN51 :N51eN46eN52

110.2e180 110.2e180 105.1e180 108.1e180 111.6e180 111.6e180

1

N4╍N10eN16╍N22eN28╍N34eN40╍N46 eN52╍N63eN68╍N73eN78╍N83eN88 N5eN11╍N17eN23╍N29eN35╍N41eN47 ╍N53eN58╍N69eN74╍N79eN84╍N89 N6╍N12eN18╍N24eN30╍N36eN42╍N48 eN54╍N59eN64╍N75eN80╍N85eN90 N7eN13╍N19eN25╍N31eN37╍N43eN49 ╍N55eN60╍N65eN70╍N81eN86╍N91 N8╍N14eN20╍N26eN32╍N38eN44╍N50 eN56╍N61eN66╍N71eN76╍N87eN92 N9eN15╍N21eN27╍N33eN39╍N45eN51 ╍N57eN62╍N67eN72╍N77eN82╍N93

2 3 4 5

distance, such as the hydrogen bond O/HeO (about 1.8 Å for O/H moiety). That is, it is the proximity nitrogen atoms at positions 10 and 4’ of the boat hexagon lead to the van der Waals interactions which are recognized with AIM analysis. From AIM analysis, It is found that there are three covalent bonds and three van der Waals interactions between two adjacent chair hexagons inserted (i.e. Chair k to Chair k þ1, 1  k  n‒2) of a molecule in (N6)n sequence (Figs. 5 and 6). Figs. 5 and 6 are threedimensional views of molecule N96 in the AIM window. Due to different angles of view, it is difficult to clearly observe the important graphic elements (covalent bonds and van der Waals forces) required. So, van der Waals forces obtained from AIM analysis are projected on 2d Fig. 1 which is the optimized structure diagram of molecule N96, and 2d Fig. 9 is obtained. Every molecule in (N6)n sequence shapes in hexagonal prism. The edges of the hexagonal prism are the connection lines interweaved with covalent bonds and van der Waals forces. This kind of lines is defined warp yarn (longitude line). The perimeter of a Chair is defined weft yarn (latitude line). According to these definitions, each molecule of (N6)n has six warp yarns and n‒1 weft yarns. The six warp yarns of the example molecule N96, i.e. (N6)16, are listed in Table 9. The average length for a warp yarn of molecule (N6)n can be estimated by the following formula: Warp length ¼ (n‒2)  (NeN bond length þ BP length)/2

(F5)

Where BP length of N40╍N46 is 2.336 Å, NeN bond length is 1.484 Å. Each warp yarn is formed by connections the covalent bond with the van der Waals, alternatively. In warp yarns, the number of covalent bonds and van der Waals is equal when n is an even number. For example, the numbers of covalent bonds and van der Waals forces are both 7 in (N6)16 which n ¼ 16. However, the number of covalent bonds is one less than that of van der Waals forces when n of the molecule (N6)n is odd. For example, the number of covalent bonds is 7 and the number of van der Waals forces is 8 in (N6)17 or N102 which n ¼ 17. Obviously, the warp yarn length calculated by (F5) is a little bit longer than it is when n is odd. So, it would make more sense to discuss the molecule structures in terms of the average warp yarn length. For the example molecule N96, i.e. (N6)16, warp yarn length of N96 ¼ (16‒2)  (2.336 þ 1.484)/2 ¼ 26.74 (Å)

6

For the largest molecule N210, i.e.(N6)35, warp yarn length of N210 ¼ (35‒2)  (2.336 þ 1.484)/2 ¼ 63.03 (Å) The nitrogen cage molecules of (N6)n are tubular. If it is taken as a nitrogen nanotube, the total length of the nanotube should be the sum of warp yarn length and the height of the two end caps. The measured height of the end cap is about 1.28 Å. So, the shortest nanotube is about 26.74 Å þ 2  1.28 Å ¼ 29.30 Å for N96, while the longest nanotube is about 63.03 Å þ 2  1.28 Å ¼ 65.59 Å for N210. The n‒1 weft yarns are the perimeters of n‒1 inserted chair hexagons of the molecule (N6)n. The atoms number on kth Chair or weft yarn are: 6k‒2, 6k‒1, 6k, 6kþ1, 6kþ2, 6kþ3 (n‒1 k  1), respectively The weft is a closed ring. As a nitrogen nanotube, the inner diameter of the tube is about 0.86 Å for all the molecules in (N6)n sequence. In Fig. 5, all the topological parameters of the example molecule N96 in AIM analysis are illustrated and the various kinds of interactions present a dense network structure. Shielding some of the parameters, the remain two are the important covalent bonds and van der Waals forces (Fig. 6 or Fig. 9) and keep its network structure. It is the existence of van der Waals that makes up the gap (No connection between atoms at position 10 and 4’ of a boat hexagon can be regarded as the openings on the prismatic surface) of the grid structure composed only of covalent bonds originally. It can be found that covalent bonds and van der Waals forces in the molecule of (N6)n interweave to form a network which tightly bind nitrogen atoms at grid points, and consequently, to make the molecules in (N6)n sequence stable. Furthermore, since a boat conformation tends to bring positions 10 and 40 closer together (atomeatom proximity), van der Waals (dispersion) interactions appear at the positions 10,40 ‒N of normal boat hexagon on a prismatic surface. The more boat hexagon on the prismatic surface has, the more number of van der Waals (dispersion) interactions, and consequently, the more stable the molecule is.

Fig. 9. The structure diagram of covalent bonds (solid lines) and van der Waals forces (dot lines) from AIM analysis of the example molecule N96 optimized at the B3LYP/cc‒pvDZ level.

Q. Guo et al. / Journal of Molecular Graphics and Modelling 96 (2020) 107508

For a pure nitrogen prismatic molecule, in principal, the surface is consisted of hexagons except the two caps. However, the prismatic surface is composed of full boat-shaped hexagon as the ends are regular triangle or quadrangle. When the ends are pentagon or hexagon, some hexagons on prismatic surface may loss boat conformation because three-coordinate nitrogen favors a pyramidal, ammonia-like environment, while the curvature of the hexagon decreases and the bonding environment of each three-coordinate atom approaches (in the infinite limit) a flat, trigonal planar arrangement. To resistant this trigonal planar arrangement, some hexagon have to convert to other form, such as chair form, etc. The decreasing number of boat conformation leads to decrease the 10,40 ‒ N proximity on the hexagons, and thus, decrease the amounts of van der Waals interactions. This is the reason why the isomers N60(D5d), and N60(D6d) are not convergent, while N60(D3d) is convergent at the level B3LYP/cc‒pvDZ in our previous work [19]. Of course, the smaller the ring on the ends, the longer the prism is. This further confirm that the van der Waals interactions increase with the length of prism, and molecules in (N6)n sequence are stable. 4. Conclusion The results of this study shows: 1) the nitrogen cages (N6)n (n ¼ 16e35) are thermodynamically stable; 2) the van der Waals forces or dispersion forces in the cage molecules play a major role in their stabilities, with a secondary stabilizing factors due to the pentagons and triangles in the cage; 3) the van der Waals forces or dispersion forces occur only on the prismatic surface; 4) it’s the van der Waals forces that makes them interweave with covalent bonds to form a dense network which tightly bind nitrogen atoms at grid points, and consequently, to make the molecules in (N6)n sequence stable; 5) the essence of the van der Waals forces or dispersion forces is the atomeatom proximity effect; 6) the boat hexagons of the atoms on the prismatic surface provide a suitable spatial configuration and orientation for the 10,40 nitrogen atoms proximity. Partial overlap of the sp3 lone pair orbitals of nitrogen atoms at positions 10,40 is the intrinsic cause of the dispersion interaction between them; 7) the more boat hexagon the prism has, the more number of van der Waals (dispersion) interactions, and consequently, the more stable the molecule is; 8) the lengths of the (N6)n (n ¼ 16e35) cages are about 2.93 nme6.56 nm, which mean that the cages may be environment-friendly beeline nanotubes. If synthesized, these nanomaterials could have broad applications.

[6]

[7]

[8]

[9] [10]

[11]

[12] [13] [14] [15]

[16]

[17]

[18]

[19]

[20]

[21]

[22]

[23]

[24] [25]

Declaration of competing interest The authors declare no conflicts of interest. Acknowledgements This work was Supported by the Foundation of Applied Basic Research Project of Sichuan Provincial Science and Technology Department (No. 2017JY0177; No. 2018JY0262).

[26]

References

[27]

[1] M.N. Glukhovtsev, H. Jiao, P.v.R. Schleyer, Besides N2, what is the most stable molecule composed only of nitrogen atoms? Inorg. Chem. 35 (1996) 7124e7133, https://doi.org/10.1021/ic9606237. [2] B.M. Gimarc, M. Zhao, Strain energies in homoatomic nitrogen clusters N4, N6, and N8, Inorg. Chem. 35 (1996) 3289e3297, https://doi.org/10.1021/ic951373h. [3] M. Tobita, R.J. Bartlett, Structure and stability of N6 isomers and their spectroscopic characteristics, J. Phys. Chem. A 105 (2001) 4107e4113, https:// doi.org/10.1021/jp003971þ. [4] Q.S. Li, D. Yong, Theoretical studies of the N6 potential energy surface, J. Phys. Chem. A 106 (2002) 9538e9542, https://doi.org/10.1021/jp0258917. [5] L.J. Wang, W.G. Xu, Q.S. Li, Stability of N8 isomers and isomerization reaction

[28]

[29]

[30]

11

of N8 (C2v) to N8 (Cs), J. Mol. Struct. 531 (2000) 135e141, https://doi.org/ 10.1016/S0166-1280(00)00438-3. A.M. Tian, F.J. Ding, L.F. Zhang, Y.M. Xie, H.F. Schaefer III, New isomers of N8 without double bonds, J. Phys. Chem. A 101 (1997) 1946e1950, https:// doi.org/10.1021/jp962878b. €tke, R.D. Harcourt, The interconversion of N12 to N8 and two T.M. Klapo equivalents of N2, J. Mol. Struct. 541 (2001) 237e242, https://doi.org/10.1016/ s0166-1280(00)00805-8. C. Chen, S.F. Shyu, Theoretical study of single-bonded nitrogen cluster-type molecules, Int. J. Quantum Chem. 73 (1999) 349e356, https://doi.org/ 10.1002/(SICI)1097-461X(1999)73:43.0.CO. D.L. Strout, Acyclic N10 fails as a high energy density, Material. J. Phys. Chem. A 106 (2002) 816e818, https://doi.org/10.1021/jp0132073. F.J. Owens, Density functional calculation of structure and stability of nitrogen clusters N-10, N-12, and N-20, J. Mol. Struct. 623 (2003) 197e201, https:// doi.org/10.1016/S0166-1280(02)00695-4. H.W. Zhou, W.X. Zheng, X. Wang, Y. Ren, N.B. Wong, Y.J. Su, A.M. Tian, A Gaussian-3 investigation on the stabilities and bonding of the nine N10 clusters, J. Mol. Struct. 732 (2005) 139e148, https://doi.org/10.1016/ j.theochem.2005.05.035. L.Y. Bruney, T.M. Bledson, D.L. Strout, What makes an N12 cage stable? Inorg. Chem. 42 (2003) 8117e8120, https://doi.org/10.1002/chin.200404005. Q.S. Li, J.F. Zhao, Theoretical study of potential energy surfaces for N12 clusters, J. Phys. Chem. A 106 (2002) 5367e5372, https://doi.org/10.1021/jp020110n. M.R. Manaa, Toward new energy-rich molecular systems: from N10 to N60, Chem. Phys. Lett. 331 (2000) 262e268, https://doi.org/10.1016/s0009-2614(00)01164-7. J. Guan, L.P. Cheng, W.G. Xu, Q.S. Li, S. Li, Z.P. Zhang, Theoretical prediction of potential energy surface for N14 cluster, J. Theor.Compu. Chem. 2 (2003) 7e14, https://doi.org/10.1142/S0219633603000331. J.D. Gu, K.X. Chen, H.L. Jiang, J.Z. Chen, R.Y. Ji, Y. Ren, A.M. Tian, N18: a computational investigation, J. Mol. Struct. 428 (1998) 183e188, https:// doi.org/10.1016/S0166-1280(97)00277-7. T.K. Ha, O. Suleimenov, M. Nguyen, A quantum chemical study of three isomers of N20 Chem, Phys. Lett. 315 (1999) 327e334, https://doi.org/10.1016/ s0009-2614(99)01271-3. D.L. Strout, Isomer stability of N24, N30, and N36 cages: cylindrical versus spherical structure, J. Phys. Chem. A 108 (2004) 2555e2558, https://doi.org/ 10.1021/jp0378889. H.W. Zhou, N.B. Wong, G. Zhou, A.M. Tian, Theoretical study on "multilayer" nitrogen cages, J. Phys. Chem. A 110 (2006) 3845e3852, https://doi.org/ 10.1021/JP056435W. H.W. Zhou, M. Beuve, F. Yang, N.-B. Wong, W.-K. Li, Theoretical investigation on the cylinder-shaped N66 cage comput, Theor. Chem. 1005 (2013) 68e74, https://doi.org/10.1016/j.comptc.2012.11.007. H.W. Zhou, N.B. Wong, G. Zhou, A.M. Tian, What makes the cylinder-shaped N72 cage stable? J. Phys. Chem. A 110 (2006) 7441e7446, https://doi.org/ 10.1021/jp062214u. H.W. Zhou, N.B. Wong, A.M. Tian, Theoretical study on the cylinder-shaped N78 cage, J. Mol. Graph. Model. 25 (2006) 578e583, https://doi.org/10.1016/ j.jmgm.2006.05.009. H.W. Zhou, N.B. Wong, Theoretical investigation on the cylinder-shaped N84 cage, Chem. Phys. Lett. 449 (2007) 272e275, https://doi.org/10.1016/ j.cplett.2007.10.076. Q. Guo, H.W. Zhou, R. Liu, Computation study on cylinder-shaped pure nitrogen molecule N90, Chem. Res. Appl. 29 (2017) 1042e1049, doi: (none). Gaussian 09, Revision C.01, M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G.A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H.P. Hratchian, A.F. Izmaylov, J. Bloino, G. Zheng, J.L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J.A. Montgomery Jr., J.E. Peralta, F. Ogliaro, M. Bearpark, J.J. Heyd, E. Brothers, K.N. Kudin, V.N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J.C. Burant, S.S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J.M. Millam, M. Klene, J.E. Knox, J.B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R.E. Stratmann, O. Yazyev, A.J. Austin, R. Cammi, C. Pomelli, J.W. Ochterski, R.L. Martin, K. Morokuma, V.G. Zakrzewski, G.A. Voth, P. Salvador, € Farkas, J.B. Foresman, J.V. Ortiz, J.J. Dannenberg, S. Dapprich, A.D. Daniels, O J. Cioslowski, D.J. Fox, Gaussian, Inc., Wallingford CT, 2009. (a) A.E. Reed, R.B. Weinstock, F. Weinhold, J. Chem. Phys. 83 (1985) 735e746, https://doi.org/10.1063/1.449486; (b) A.E. Reed, L.A. Curtiss, F. Weinhold, Chem. Rev. 88 (1988) 899e926, https://doi.org/10.1021/CR00088A005. (a) F. BieglereKünig, J. Schünbohm, R. Derdau, D. Bayles, R.F.W. Bader, AIM 2000, Version 2.0, McMaster University, Hamilton, 2002; (b) R.F.W. Bader, Atoms in Molecules, A Quantum Theory Oxford University Press, Oxford, 1994, pp. 275e426. J.M. Molina, J.A. Dobado, M.C. Daza, J.L. Villaveces, Structure and bonding of weak hydrogen peroxide complexes, J. Mol. Struct. 580 (2002) 117e126, https://doi.org/10.1016/s0166-1280(01)00602-9. P.L.A. Popelier, Fast algorithm to compute atomic charges based on the topology of the electron density, Theo. Chem. Acc. 105 (2001) 393e399, https:// doi.org/10.1007/s002140000224. U. Koch, P.L.A. Popelier, Characterization of C-H-O hydrogen-bonds on the basis of the charge-density, J. Phys. Chem. 99 (1995) 9747e9754, https:// doi.org/10.1021/j100024a016.