Topological and Quantum Stability of Low-Dimensional Crystalline Lattices with Multiple Nonequivalent Sublattices

The terms of topological and quantum stabilities of low-dimensional crystalline carbon lattices with multiple non-equivalent sublattices are coined using theoretical analysis, multilevel simulations, and available experimental structural data. It is demonstrated that complex low-dimensional lattices are prone to periodicity breakdown caused by structural deformations generated by PBC approach. To impose PBC limitations for low-dimensional lattices, the Topology Conservation Theorem (TCT) is formulated and proved. It is shown that the lack of perfect filling of 2D crystalline space units may cause formation of i) Structure waves; ii) Nanotubes or rolls; iii) Saddle structures; iv) Aperiodic atomic clusters; v) Stabilization of 2D lattices by aromatic resonance, correlation effects, or van-der-Waals interactions. The effect of quantum instability of infinite structural waves is studied using quasiparticle approach. It is found that both perfect finite-sized, or stabilized structural waves can exist and can be synthesized. It is shown that for low-dimensional lattices prone to breakdown of translation invariance (TI), complete active space of normal coordinates cannot be reduced to a subspace of TI normal coordinates. As a result, constrained TI subspace structural minimization may artificially return a regular point at the potential energy surface as either a global/local minimum/maximum. It is proved that phonon dispersion cannot be used as solid and final proof of their stability. It is shown that ab initio molecular dynamics (MD) PBC Nos\'e-Hoover thermostat algorithm constrains the linear dimensions of the periodic slabs in MD box preventing their thermostated equilibration. Based on rigorous TCT analysis, a flowchart algorithm for structural analysis of low-dimensional crystals is proposed and proved to be a powerful tool for theoretical design of advanced complex nanomaterials.


Introduction
Existence of low-dimensional lattices is limited by fundamental thermodynamic instability of one-(1D) and two-dimensional (2D) crystals at finite temperatures (Landau-Peierls [1] and Mermin-Wagner [2] theorems). Both theorems state that for extended 1D and 2D crystals the crystallographic long-order cannot be maintained because of logarithmic unconvergency of mean square displacements of lattice atoms even at T = 0 (1D case) or at the lower limit of integration (2D case). According to the theorems, long-wavelength fluctuations destroy the long-range order of low-dimensional crystals. 2D lattices embedded in a 3D space have a tendency to be crumpled [3] by temperature fluctuations. The anharmonic coupling between bending and stretching modes leads to suppression of temperature fluctuations finally allowing 2D crystals to be exist exhibiting strong height fluctuations [3][4][5].
Theoretical predictions of the atomic and electronic structure and properties of key nanoclusters have played a leading role in the development of a new revolutionary class of low-dimensional (0D, 1D and 2D) materials. Prediction of fullerenes C 60 [6,7] and C 20 [7] was made back in 1970 and 1973, respectively, using simple structural consideration of aromatic systems [6] and Hückel electronic structure calculations [7]. Experimental discovery of C 60 buckminsterfullerene in 1985 [8] excellently confirmed all predicted features of atomic and electronic structures of the fullerenes.
The layered crystalline structure of graphite was well known long before experimental exfoliation of graphene [9], so one can consider the first theoretical investigation of graphene electronic structure made by tight-binding model [10,11], linear combination of atomic orbitals (LCAOs) tight-binding [12] and local density approximation (LDA) [13] within periodic boundary conditions (PBCs) as pioneering studies of band structures of 2D crystalline lattices. Later, theoretical results were confirmed with high accuracy by the ARPES technique [14].
Just after the famous discovery of carbon nanotubes [15], the electronic structure of the tubulenes was successfully independently predicted using tight-binding [16,17] and DFT [18,19] approximations within the PBC approach. Later, theoretical predictions were confirmed by scanning tunneling spectroscopy studies [20].
Inspired by discovery of carbon nanotubes, the atomic and electronic structure of h-BN nanotubes were predicted and studied theoretically back in February [21] and November [22] of 1994 using tight-binding and LDA calculations within PBC approximation. Experimental discovery of the tubes followed shortly in 1995 [23,24] was a great experimental achievement based on brilliant theoretical prediction.
Another fascinating prediction of structure and properties of diamanes, the thinnest diamond films, which consist of two 111 cubic or 001 hexagonal diamond planes, respectively, was also made based on density functional (DFT) plane wave (PW) PBC calculations [25]. Just recently, following predicted reaction pathway fluorinated single-layer diamond film was synthesized by fluorination of bilayer graphene [26].
Biphenylene-based novel 2D carbon π-conjugated material [27] composed of a combination of four-, six-and eight-membered carbon fragments, as well as nanoribbons and nanotubes on its base was proposed and studied using ab initio DFT simulations. The low-dimensional atomic lattices were predicted to demonstrate promising hole/electron mobility values typical for ambipolar organic semiconductors with band-gap being significantly smaller than that for the 1D polymer ribbons containing the same number of biphenylene units. The bottom-up approach [28] was used for growth of 2D biphenylene sp 2 -hybridized carbon network through an interpolymer lateral dehydrofluorination fusion of benzenoid polyphenylene chains on top of Au (111) surface with consequent C-C coupling to form a biphenylene regular lattice. It was experimentally characterized by scanning probe methods, which revealed its metallic nature. The synthesis was accompanied [28] by PBC DFT electronic structure calculations, dI/dV differential conductance spectra and theoretical simulation of biphenylene network nc-AFM images using 9 × 1 supercell. Comparison of experimental structural data and PBC DFT calculations demonstrated perfect coincidence of experimental and theoretical results. One should mention theoretical prediction [29] and followed experimentally discovery [30] of 2D borophenes as well.
An existence of standing waves in graded 1D InAs/GaAs core-shell nanowires under temperature gradient was predicted using nonequilibrium molecular dynamics simulations [31], and it was found that the standing wave greatly enhance nanowire thermal rectification effect. Using non-equilibrium molecular dynamics structural waves were predicted as well in (0001) oriented wurtzite InAs nanowire and zinc-blend (110)-oriented Si nanowire [32]. In particular, it was shown that the length of the nanowires determines whether or not a standing wave can be formed and the nanowire cross-sections can influence the natural frequencies along the lateral directions. Using nonequilibrium molecular dynamics simulations [33] standing waves also were predicted for non-equilibrium steady-state for graphene and h-BN superlattices. A detailed analysis of the phonon spectra shows that this anomalous thermal conductivity behavior is a result of strong phonon wave interference.
As mentioned above, most theoretical predictions were made using electronic structure calculations based on the PBC approach. The PBC approach is a powerful and widely used tool in solid-state physics and material science to study atomic and electronic structure and spectroscopic characteristics of a vast variety of crystalline lattices, including superlattices. In fact, the PBC approximation can be considered as mandatory symmetry operation, which generates linear 1D, planar 2D and full space 3D crystalline lattices. Since PBC is a kind of symmetry restriction, one can speculate about its improper application for calculations of low-dimensional crystalline lattices with non-zero curvature, which may result in artificial linear crystalline lattices.
Phonon spectra calculation is a well-known, popular and mandatory test of structural stability of the lattices studied within the PBC approach (see, for example, [34]), where the absence of imaginary modes (all force constant matrix elements are non-negative) is considered as a solid and final proof of stability of proposed crystals. It is presumed that in crystals the atoms move around their equilibrium positions following harmonic law. The harmonic approximation works well for 3D crystalline lattices in the vicinity of structural global minima, while obviously it is not the case for 2D/1D materials when one can easily consider a number of elastic modes with linear or even constant dependence of the lattice total energy upon the small bending or twisting deformations perpendicular to the main crystalline plain/axis due to negligibly small corresponding force constants.
In recent years a vast number of low-dimensional nanomaterials were predicted with promising properties based on theoretical PBC calculations (among others see, for example, [35,36]). Penta-graphene [35] is likely one of the most famous members of the penta-family made of only pentagons. It is characterized by a three-layer thick structure and resembles Cairo pentagonal tilling with central sp 3 carbon sublattice and two perpendicular symmetrically nonequivalent sp 2 carbon sublattices located in the top and bottom layers. The constituting elements can vary and may include Si, N, and other either metallic or non-metallic ones. In contrast with graphene, which is, in fact, one-atom-thick material, penta-graphene has non-zero thickness (1.2 Å). Penta-graphene was theoretically proposed [35] to be thermodynamically and mechanically stable carbon allotrope, which stability was proved as well by phonon spectrum calculations (no imaginary modes) [35]. However, further theoretical investigations showed that it is mechanically [37], thermodynamically [38], symmetrically and topologically [39,40], energetically and kinetically [41] unstable. In particular, it was shown [39], that mutually perpendicular top and bottom sublattices of penta-graphene create uncompensated mechanical stress, which causes strong topological instability and bending of the 2D crystalline lattice. Though some unstable materials can be stabilized by a suitable substrate but this is not the case of penta-graphene because of large internal stress [39] which cannot be compensated by van-der-Waals interactions by any kind of supports.
Another perfectly planar 2D carbon allotrope, so-called 'phagraphene' was predicted to be stable in reference [36] with rectangular unit cell parameters a 1 = 8.0946 Å and a 2 = 6.6518 Å. Its crystalline lattice consists of a regular combination of 5-, 6-and 7-member carbon rings with each single pentagon-heptagon pair surrounded by 8 carbon hexagons. The stability of phagraphene was proved by calculations of phonon spectra [36,42] (no imaginary modes) using PHONON [43][44][45] and GULP [46,47] codes, respectively. Later, using tight-binding and density functional theory it was shown [48] that phagraphene has at least two possible configurations, namely planar and wavy one, which are nearly degenerate in energy.
Different disclinations (see, for example, [49][50][51][52]) and regular 2D sp 2 carbon lattices composed of a combination of 5-, 6-, and 7-member carbon rings were proposed long before phagraphene [36] (some of the lattices will be considered below). Regular planar 2D Stone-Wales (SW) lattice [49,53,54] (pentaheptite or R 5,7 Haeckelite) composed of only carbon pentagons and heptagons was proposed and the atomic and electronic structure was studied using tight-binding and LDA PW PBC approaches. The R 5,7 can be generated by the rotation of each single carbon-carbon bond by 90 • angle in the center of four fused carbon hexagons of graphene lattice. Significantly more complex planar 2D regular lattices composed of 5-, 6-, and 7-member rings, so-called Haeckelites (namely hexagonal H 5,6,7 , and oblique O 5,6,7 lattices) were proposed and studied using the tight-binding approach in reference [49]. Among three types of Haeckelites, the hexagonal one, H 5,6,7 , is of special interest because it is composed of linear 5/6/7 cores arranged in regular tilling of hexagonal symmetry. Other regular planar one-atom-thick lattices of intercalated compounds C 6 Ca and C 6 Yb composed by 5-, 6-, 7-, and 8-member carbon rings under high pressure (>18 GPa) were studied [55][56][57] using the DFT PW PBC approach and it was found that the most stable Cmmm lattice (which was obtained by rotation of two carbon-carbon bonds) is composed of a combination of five-and eightmember rings and the next in energy Pmma lattice is composed by lines of 5/7 cores (which is in fact are lines of fused SW defects) separated by lines of carbon hexagons of one hexagon width. Planar and out-of-plain [52,[58][59][60][61][62]] sp 2 conjugated carbon lattices with SW defects were also considered as well.
One can easily find dozens of predicted unusual/fascinating/fantastic one-, two-and three-atom thick 2D crystalline lattices composed by 3-, 4-, 5-, 6-, 7-, 8-, and even 9-member rings of different light and heavy elements like boron, carbon, nitrogen, silicon, different metals etc somewhere in the literature (in this work the authors are not intended to make the references on disputed publications without any obvious urgent need, nor all of the mentioned hypothetical materials will be critically considered in this study). Mostly the predictions were made using PBC calculations at different levels of theory with calculated phonon spectra as solid and final proof of the stability of the lattices.
Using atomistic simulations [63] it was found that in pristine narrow graphene ribbons, disruption of the aromatic bond network results in depopulation of covalent orbitals and tends to elongate the edge, with an effective force of f e ∼ 2 eV Å −1 . In the case of narrow ribbons, this effect favors structural spontaneous twisting, resulting in the parallel edges forming a double helix with a pitch of about 15-20 lattice parameters. For wide graphene ribbons the twist abruptly vanishes and instead the corrugation localizes near the edges. Using elastic plate theory [64], it was found, that the edge stresses at both zigzag and armchair edges introduce intrinsic ripples in freestanding graphene sheets even in the absence of any thermal effects. It causes out-of-plane warping to attain several degenerate mode shapes. The scaling laws for the amplitude and penetration depth of edge ripples as a function of wavelength were identified. It was found that the edge stresses can lead to twisting and rolling of the nanoribbons.
In this work, a theoretical model to consider mechanical tension in low-dimensional lattices was proposed and it was found that multiple nonequivalent sublattices may generate internal structural stress caused by uncompensated torques. Based on symmetry and structural analysis, a topology conservation theorem (TCT) was formulated and proved for 1D and 2D crystalline lattices. According to TCT, to avoid the uncompensated mechanical stress perpendicular to the main plane or axis of a low-dimensional lattice, the constituting structural fragments must perfectly fit the low-dimensional space without implying any external forces. Due to the leading contribution of the stretching force constants to total energy, any small structural mismatch leads to bending of low-dimensional crystalline lattices in the perpendicular direction to the main plain or axis to eliminate the resulting mechanical stress. It was shown that PBCused to introduce linear translation invariance (TI) symmetry restrictions may cause artificial structural stabilization of linear 1D and planar 2D lattices which are prone for structural deformations. This may lead to multiplication of the translation periods with formation of superperiodic structural waves or even breakdown of the lattice TIin the forms of screwing, bending, or rolling. In the worst case, the PBC may artificially stabilize low-dimensional lattices which are otherwise structurally unstable with breakdown of fundamental short-order structural parameters like bond symmetry, coordination numbers and nature of chemical bonding itself like hybridization character. It was found that in many cases the energy of artificial structural stabilization surpasses the energy of van-der-Waals interactions, which makes it impossible to stabilize the low-dimensional nanoclusters on any kind of substrate and consequently synthesize them. Cluster optimizations of different sizes of finite flakes at various levels of theory result either in bent lattices with multiplication of translation vectors, breakdown of periodicity, structure waves of variable or constant wave lengths, 1D nanotubes or rolls, or even barrier-free transformations to separate graphene-like flakes or aperiodic irregular carbon clusters. Hessian calculations of finite clusters constituted by 5-, 6-and 7-member rings revealed imaginary vibrational frequencies of all considered planar conformers, which unequivocally prove their transition state nature. PBC optimizations of 1D nanoribbons of one-atom-thick wavy carbon allotropes confirm significant energy gain of wave-like structures in comparison with planar conformers. It was shown that structure of standing waves of finite-length flakes can be described by homogeneous D'Alamber second-order differential equations with predetermined boundary conditions. The structure of infinite free-standing lattice waves can be described by inhomogeneous D'Alamber second-order differential equations, for which the solutions can be written as infinite decomposition over planar waves with non-zero external forces acting on the lattice atoms. In contrast to infinite free-standing wave lattices, the finite lattice waves were considered as quantum-stable structures due to perfectly localized positions of the atoms. The solutions of inhomogeneous equations for free-standing infinite lattices were interpreted in terms of quantum instability due to complete delocalization of atomic positions of lattice constituting atoms. It was shown that artificial stabilization and lack of imaginary modes in phonon dispersions of planar conformers prone for breakdown of periodicity may be caused by constrained structural minimization in a subspace of TI normal coordinates (TI SS NC) of true Hilbert-space potential energy surfaces far away from any kind of extremum points. It consequently leads to improper localization of a regular point in complete active space of normal coordinates (CAS NC) as a potential energy surface extremum with consequent erroneous phonon dispersion law calculated at constrained extremmum. It was concluded that the absence of imaginary modes in phonon dispersion laws cannot be considered as a solid and final proof of either structural stability or metastability of low-dimensional crystalline lattices with multiple nonequivalent sublattices. Almost perfect coincidence of the cluster optimizations at various levels of theory combined with PBC PW DFT consideration completely excludes the electronic factor in artificial structural stabilization and directly indicates the symmetry nature of the effect. It was found that in some special cases the structural stress caused by departure of the lattice from mandatory TCT requirements can be almost completely compensated by aromatic resonance and correlation effects. It was shown that ab initio molecular dynamics (MD) PBC Nos´e-Hoover thermostat algorithm constrains the linear dimensions of the periodic slabs in MD box preventing their thermostated equilibration. Based on TCT analysis, a flowchart algorithm to study structure and properties of low-dimensional crystals was proposed and proved to be a powerful tool for theoretical design of advanced complex nanomaterials. Proposed theoretical results perfectly correspond to available structural experimental data.

Structural models and computational methods
To study artificial structural stabilization of low-dimensional lattices caused by translation symmetry restrictions imposed by PBC, one-, two-and three-atom thick low-dimensional carbon sp 2 and sp 3 allotropes with multiple nonequivalent sublattices were designed by regular combination of rectangular, five-, six-, seven-, and eight-member structural fragments. Along with a set of complex low-dimensional carbon crystals with multiple sublattices, a number of finite graphene flakes (figure 1(a)) were considered as a reference to check the accuracy of atomic and electronic structure calculations. It is necessary to note that in this study a set of low-dimensional carbon sp 2 /sp 3 lattices were chosen for simplicity since it is well known that symmetry and atomic structure of the lattices are mostly determined by hybridization nature of constituting carbon atoms which allows one to easily separate symmetrical and electronic structural contributions.
The SW lattice [49,53,54] which can be obtained from graphene lattice through C-C bonds rotation caused by SW transformation [53] is presented in figure 1(b). The SW transformation converts all hexagons of graphene lattice into pentagon/heptagon pairs, effectively eliminating all graphene hexagons. In R 5,7 Haeckelite all pentagons and heptagons (and 5/7 pairs as well) share common C-C edges.
The 'phagraphene' lattice [36] is presented in figure 1(c). It was designed by a combination of planar five-, six-and seven-member carbon rings. The main structural feature of phagraphene is pairs of carbon penta-and heptagons which share one common carbon-carbon bond (5/7 fragment), surrounded by 6 hexagons. The neighboring 5/7 pairs contact each other by the vertices of pentagons and heptagons through one C-C bond.
Another phagraphene-like lattice can be designed by insertion of one hexagon between pentagons and heptagons effectively separating the 5 and 7 fragments (figure 1(d)). Let us introduce p(pentagon) h(heptagon)C(carbon), phC(n, m), notations to denote 2D lattices composed of 5-, 6-, and 7 carbon fragments, where n index denotes the number of hexagons, which separate 5 and 7 rings, and m index indicates the number of hexagons which separate the 5/7 pairs. In these notations, the phagraphene lattice can be denoted as phC(0,1), the SW lattice, R 5,7 (no hexagons in the lattice), must be indexed as phC(0,0), the next one, with one hexagon separating pentagons and heptagons and without C-C bonds separating 5/6/7 cores can be denoted as phC (1,0), and the last considered lattice with completely separated pentagons and heptagons by the hexagons (each single pentagon surrounded by 5 hexagons and each single heptagon surrounded by 7 hexagons) is denoted as phC(1,1) (figure 1(e)). All four phC(0,1), phC(0,0), phC (1,0), and phC (1,1) lattices are investigated in details in this study. One can simply imagine an infinite number of regular phC(n,m) lattices with arbitrary large n and m indexes (see, for example [59] for some of such structures).
The one-atom-thick 2D biphenylene lattice [27,28] is composed by fused biphenylene fragments and is a combination of four-, six-and eight-membered carbon rings (figure 1(f)). Each single elongated hexagon in 2D biphenylene lattice is surrounded by 4 elongated planar octagons and 2 four-member carbon rings rather close to regular squares. Each single carbon square is surrounded by 2 hexagons and 2 octagons and the octagons are surrounded by 2 other octagons, 2 squares and 4 hexagons. It is necessary to note that despite of very specific electronic effects which will be described below, the lattice can be considered as a planar low-dimensional crystal.
Three-atom-thick mixed sp 2 /sp 3 penta-graphene (figure 1(g), [35]) which consists of a combination of only pentagonal fragments was tested as well for artificial structural stabilization caused by PBC translation symmetry restrictions. Penta-graphene was the first among many others three-atom-thick (one can simply find dozens of publications devoted to theoretical predictions of penta-like 2D materials) theoretically 'discovered' 2D crystalline lattices constituted only by pentagonal fragments, which stability was 'confirmed' by phonon spectra calculations in the framework of translation symmetry restrictions imposed by linear PBC approximation.
The last regular hypothetical 2D carbon allotrope, ladderane [65], consists of sp 3 atoms arranged in a two-atom thick buckled rectangular lattice (figure 1(h)). Ladderane can be characterized by two different types of four-coordinated C-C bonds. The first type consists of two distorted single carbon-carbon bonds with perfect 180 • angle between them which follow linear pattern forming perfect atomic lines which belong to a single atomic plane (two atomic planes in total). Two other single C-C bonds keep close to perfect sp 3 bond angle between them bonding atomic lines from neighboring atomic planes. In total, each single carbon atom of ladderane lattice has two single C-C bonds with the partners in the atomic line and two bonds which connect the atomic planes.
The finite flakes and nanoribbons cut out from parent 2D crystals with multiple sublattices were considered using the cluster approach. First, to elucidate the role of symmetry restrictions and exclude the electronic factors in determination of the main features of the lattice structure, the model MM+ potential (a development of the MM2 method [66]) was used to simulate the atomic structure of finite extended clusters. Semiempirical PM3 [67] and hybrid DFT B3LYP [68,69] coupled with Pople split-valence 6-31G * basis set [70,71] lattice as implemented in the GAMESS package [72] were used for structural optimization of finite carbon clusters. For some finite clusters of phC(1,0) phC(0,0) and phC (1,1) lattices, combination of 3-21G and 6-31G * basis sets was used as well. Model potential and semiempirical PM3 approaches were used to calculate the atomic structures of extended (up to 1000 carbon atoms) nanoclusters and to reveal the role of electronic factors and accuracy of structure determination as references to ab initio DFT B3LYP approach. All cluster electronic structure calculations were performed without taking into account Grimme corrections [73] for van-der-Waals interaction since disperse interactions is an extence characteristic which makes impossible to comare the results of electronic structure calculations for the clusters of different sizes.
The atomic and electronic structure and vibration spectra of a number of graphene flakes (C 42 H 18 ,  C 114 H 30 , C 222 H 42 , C 366 H 54 , C 546 H 66 , and C 762 H 78 ) of hexagonal symmetry with armchair edges (figure 1(a)) were calculated to test reliability of Hessian calculations of extended low-dimensional carbon clusters at PM3, ab initio B3LYP/3-21G and B3LYP/6-31G * levels of theory. It was shown that Hessian calculations of all small flakes (up to PM3 C 222 H 42 ) returned only real vibration frequencies starting from 0.01 cm −1 as the smallest ones. Depending on the accuracy of structural optimization (which should be 10 −5 kcal/mol/Å or even better for Hessians making it difficult to achieve), for extended flakes (C 366 H 54 , C 546 H 66 , and C 762 H 78 ), Hessian PM3 calculations may return one or two imaginary modes of ∼2i cm −1 , that may be considered as practical accuracy of vibrational mode calculations.
The PBC calculations were performed using Vienna ab initio simulation package (VASP) [74,75] coupled with projector augmented wave (PAW) [76] basis set and GGA PBE functional [77] with wavefunction cutoff energy equal to 400 eV. The first Brillouin-zone for periodic calculations was sampled depending on structural features according to the Monkhorst-Pack scheme [78]. The convergence tolerances of forces and electronic minimizations were 10 −3 eV Å −1 and 10 −5 eV, respectively. All periodic images were separated by a vacuum distance at least of 15 Å to avoid artificial interactions of the lattice images along z-direction. Several finite clusters were calculated at Γ-point of the Brillouin zone. The visualization for electronic and structural analysis software (VESTA) [79] was used for visualization of the results. For selected 2D and 1D crystalline lattices, the PHONOPY code [34] was used to calculate phonon spectra. For all PBC calculations weak disperse interactions were taken into account using Grimme empirical corrections [73].
In contrast to PBC optimization procedure, during which both atomic coordinates and translation vectors are minimized, all thermostat algorithms conserve topology of simulated ensembles. In particular, NVT algorithm keeps the volume of MD box constant, keeping constant translation vectors, and recalculating atomic coordinates at each step of MD simulation.Effectively, it freezes the MD box dimensions preventing the dimension changes of 2D and 1D objects confined in the box. MD simulations of 0D objects like molecules and clusters (it could be nanocrystals and nanoflakes, for example), do not have such disadvantage since finite atomic clusters are not strictly confined by fixed MD box dimensions, which allow them to change freely their dimensions and orientation.

1D h-BN zig-zag narrow nanoribbon case
A hypothetical 1D h-BN zig-zag narrow nanoribbon (h-BN ZNR) of one B 3 N 3 hexagonal fragment width (figure 2) may be considered as the simplest case of a low-dimensional lattice with multiple non-equivalent sublattices. Without loss of generality one can introduce a 1 , a 2 , and b 1 , b 2 nonequivalent sublattices for N and B basis atoms, respectively, associated with symmetrically nonequivalent a a1 , a a2 , a b1 , a b2 translation vectors oriented along X direction. In particular, both external atomic rows have coordination numbers equal to 2 with symmetrically different environment, and different neighborhood, namely a 1 N atoms have 2 b 2 B neighbors while b 1 B atoms have 2 a 2 N neighbors. The coordination numbers of a 2 N and b 2 B atoms are equal to 3 with 2 b 1 and 1 b 2 boron atoms and 2 a 1 and 1 a 2 nitrogen neighbors, respectively. For simplicity all interatomic bond lengths can be considered equal to R NB . Since the basis atoms, the boundary conditions, coordination numbers and types of environment of a 1 /a 2 and b 1 /b 2 N and B sublattices are nonequivalent, one can write: a a1 = a a2 = a b1 = a b2 . Figure 2. A fragment of 1D h-BN zig-zag narrow nanoribbon of one B 3 N 3 hexagonal fragment width, oriented along X-axis. a 1 , a 2 , b 1 , and b 2 are nonequivalent sublattices filled by N (a 1 , a 2 , in blue) and B (b 1 , b 2 , in green) atoms, respectively. R NB is the length of interatomic bonds, a a1 , a a2 , a b1 , a b2 are translation vectors associated with a 1 , a 2 , b 1 , b 2 sublattices, respectively. α is the angle between X direction and r, where r is an atomic radius vector.
The N a 1 and B b 1 sublattices with coordination numbers equal to 2 possess different force constants Q 1 , Q 2 of vibrations along X direction because of different symmetry and environment of the edge N and B atoms even a 2 -b 1 -a 2 and b 2 -a 1 -b 2 angles and R NB length of all N-B bonds keep initial parameters of perfect 2D h-BN. The force constant Q along the X direction and R NB length of B-N bonds for a 2 and b 2 sublattices with coordination numbers equal to 3 keep the initial Q i value of perfect 2D h-BN lattice. For perfect hexagonal lattice |a a1 | = |a b1 | = √ 3R NB . So, for a a2 and a b2 translation vectors and the force constants the following relationships can be written: In harmonic approximation and taking into account symmetry restrictions, for the stress energies E i and forces acting on a i and b i sublattices one can write: Symmetrical non-equivalency of the forces acting on the boundary B and N atoms (F a1 = F b1 ) creates mechanical stress and structural curvature of ultranarrow h-BN ZNR and other similar zig-zag heteroatomic nanoribbons [82].
For perfect hexagonal lattice the torques τ i caused by uncompensated forces acting on nonequivalent sublattices of h-BN ZNR in respect to the motionless center of the mass can be written as and a 2 , b 2 associated torques mutually compensate each other. For a 1 , b 1 sublattices of the opposite N and B atoms the α angles are equal to −π/2 and π/2, respectively, (figure 2), and The τ b1 and τ a1 are oriented in opposite directions. The total torque acting on the nanostructure is the sum of the torques associated with the sublattices. Since F a1 = F b1 and the torques do not compensate each other. Since both nonequivalent torques associated with a 1 and b 1 sublattices are perpendicular to the plane of h-BN ZNR lattice and oriented in opposite directions, they lead to displacements of the atoms with the formation of a cone fragment of one B 3 N 3 width and breaking down the linear translation symmetry [39,82]. Detailed derivation of the formulas are presented in Supplemetary Data SI2 section.

1D graphene zig-zag narrow nanoribbon case
The model can be easily extended for the case of 1D zig-zag graphene narrow nanoribbon, which has equivalent edges with equivalent force constants Q 1 , Q 2 , so the associated torques perfectly compensate each other and the nanoribbon keeps its linear translation symmetry.

2D graphene and h-BN cases
Both 2D graphene and h-BN crystalline lattices have two basis atoms in hexagonal unit cells associated with 2 translation vectors a 1 and a 2 with absolute values |a 1 | = |a 12 | = √ 3R, where R is the length of either C-C or B-N chemical bonds. The angle between translation vectors, α = π/3 or α = 2π/3 can be chosen in two different ways depending on the choice of the unit cell. For the sake of simplicity let us consider the case of the first graphene sublattice (figure 3) with α 1 = π/3 and the force constant Q.
Because of the symmetry of the lattices, the partial energy components acting on atoms 1, 2, 3 are equal to each other: which can be satisfied only if all components are equal to 0, so, for pristine hexagonal lattices the X and Y components of the forces mutually compensate each other and the total stress of the lattice is equal to 0. For the torques acting on atoms 1 and 2 in respect to the origin of the coordinate system one can write: and the torques mutually compensate each other as well, keeping TI and planar 2D structure. Detailed derivation of the formulas are presented in SI2 section.

Inclusion of pentagons and heptagons into hexagonal carbon lattice
According to Theorema egregium, geometrically, an introduction of one pentagonal fragment in perfect hexagonal lattice leads to formation of positive 30 • Gaussian curvature [83,84], whereas introduction of heptagon leads to negative 30 • Gaussian curvature [85][86][87][88][89][90][91][92]. It is necessary to note that both planar and curved conformers of all low-dimensional carbon clusters considered in this paragraph were optimized by MM+, PM3 and ab initio B3LYP/6-31G * methods without applying of any symmetry restrictions. The C 20 H 10 corannulene or [5]circulene (figure 4(a)) [93] and C 28 H 14 [7]circulene [94,95] (figure 4(b)) are the well-known examples of the molecules with distinctively pronounced convex and saddle shapes of C 5v and C s symmetries, respectively, in their ground states. For both of these molecules, two equivalent curved global minima are separated by planar D 5h and D 7h transition states, respectively. Fusion of one pentagon results in formation of a perfect planar C 10 cluster (5/7 core) of C 2v symmetry. Consequent adding of 8 C 6 rings around 5/7 core leads to formation of C 32 cluster (5/7 core surrounded by 8 hexagonal rings, 5/7/6 flake, figure 4(c)) with the C 2v symmetry of planar conformer. One could speculate that fusion of +30 • and −30 • curved clusters should lead to cancellation of the curvature of the resulting flake with restoration of perfectly planar sp 2 carbon lattice. Initial C 10 and C 32 flakes were designed keeping the structure and symmetry of 5/7 core and the lengths of all C-C bonds equal to 1.4 Å. In particular, the pentagonal fragment of the 5/7 core is characterized by internal 108.00 • angles between the lattice nodes, as well as the heptagonal fragment keeps original 128.57 • . The external angles between the nodes at 5-7 boundary are equal to 123 • . One can easily discover significant departures of the angles in the hexagons which surround the 5/7 core (figure 4(c)) from perfect 120 • to 104 • -132 • . The angle between the lines connected external 3 nodes of both sides of the C 32 cluster is equal to 21 • . Significant departure of hexagon internal angles in the 5/7/6 flake from the parent hexagonal lattice angle of 120 • makes impossible perfect filling of the hexagonal lattice by introducing the 5/7 defect. Optimization of C 32 H 14 5/7/6 flake at model potential and B3LYP/6-31G * levels of theory keeps its perfect planar structure of C2v symmetry.
Semiempirical PM3 optimization of 5/7/6 C 32 H 14 flake revealed planar C 2v as well as two equivalent curved C s conformers (figure 4(c)), which display distinctive uncompensated curvature of complex nature caused by mechanical stresses introduced by 5/7 core. At the PM3 level of theory, curved conformers are symmetrically equivalent global minima without imaginary frequencies in vibrational spectrum and the smallest real frequency 28.43of cm −1 ). Planar conformer is a transition state with the barrier 0.418 kcal mol, which is much smaller than the accuracy of the PM3 method [96]) between the minima with one imaginary vibration frequency of 39.78i cm −1 and first real frequency of 22.00 cm −1 . Comparing PM3 relative energies, mechanical stress of 5/7/6 C 32 H 14 cluster caused by interaction of 5-and 7-fragments competes in energy with aromatic resonance effect, making both conformers almost equal in energy. Here and after the PM3 method is used to reveal the role of electronic factors taking into account at different levels of rigidity in determination of fine details of atomic structure of low-dimensional lattices prone to deformations caused by internal mechanical stress. At restricted R-B3LYP/6-31G * level of theory the 5/7/6 C 32 H 14 cluster retains perfectly planar C 2v symmetry which might indicate that the π-system resonance energy exceeds the deformation energy conserving perfectly planar cluster geometry. Comparison of closed-shell R-B3LYP/6-31G * and opened shell U-B3LYP/6-31G * singlets and triplet states returned closed-shell singlet as the ground state with 35.317 kcal/mol singlet-triplet splitting.
Introduction of 5/7 core in extended graphene flake like C 522 H 56 cluster with 8 C 6 rows around the central 5/7 C 10 fragment (figure 4(d)) leads to departure the symmetry from parent perfect D 6h of graphene flake to C s point group caused by insertion of one additional carbon hexagon to both left and right bottom edges of the flake with formation of regions with uncompensated positive and negative curvatures (the number of C 6 layers can be either increased or decreased without loss of generality). At the model potential level of theory, the perfectly planar transition state has relative energy E Rel planar = 51.343 kcal mol in respect to two equivalent curved C 522 H 56 global minima with relative energy E Rel curved = 0.0 kcal mol.

Structural units and fragments of phC(0,1) lattice with pair of antiparallel 5/7 cores
The unit cell of phagraphene [36] (or in terms introduced above, the phC(0,1) lattice) consists of two antiparallel 5/7 fragments separated by 3 carbon hexagons (the smallest phC(0,1) cluster constituted by two 5/7 pairs separated by 3 hexagonal rings with additional 3 + 3 hexagons from left and right, C 40 , is presented in figure 5(a)). The unrelaxed C 50 planar fragment of phagraphene (figure 5(b)) with additional surrounding 10 atoms was designed keeping all C-C bonds equal to 1.4 Å and bond angles of the 5/7 core (figure 4) of 108 • and 128 • for pentagons and heptagons, respectively. The hexagon between 5/7 partners (just in the center of figure 5(b)) which separates 5/7 partners keeps the same angles as the most distorted hexagon in C 32 cluster of single 5/7 core (figure 4(c)), namely 105 • , 123 • and 132 • which strongly depart from 120 • angles of perfect hexagons. In spite of strong distortions of hexagons, all external zigzag edges of the unit cell are perfectly parallel to each other due to squeezing of the central hexagon between 5/7 cores and stretching of the external sides of its two adjusting hexagons. Comparison of unrelaxed phagraphene C 50 fragment with the perfect hexagonal lattice (figure 5(b)) reveals significant multiple structural departures up to 0.919 Å of the 5/7 cores along both zigzag and armchair directions of graphene lattice in addition to the distortions of the hexagons which constitute the phagraphene fragment. In particular, at the equator, the unrelaxed phagraphene unit cell is 1.428 Å wider than the parent hexagonal lattice. Optimization of C 40 H 16 phagraphene fragment (figure 5(a)) at model MM+ potential level of theory returned planar cluster of C 2h symmetry. At the PM3 level of theory two symmetrically equivalent curved global minima conformers of C 2 symmetry with relative energy E Rel curved = 0.0 kcal mol and planar transition state of C 2h symmetry (two imaginary frequencies of 71.31i and 69.97i cm −1 and the first real frequency of 51.07 cm −1 ) with relative energy E Rel planar = 12.758 kcal mol were located. The ab initio R-B3LYP/6-31G * also returned two curved C 2 equivalent global minima (E Rel curved = 0.0 kcal mol) and one planar transition state of C 2h symmetry with relative energy E Rel planar = 6.233 kcal mol. Improper torsion angles along and perpendicular to 5/7 cores for PM3 curved global C 40  To examine the perfect environment of 5/7 cores of phagraphene, the C 50 H 18 cluster was considered (figure SI1-1), in which two central 5/7 cores are separated by three C 6 rings and surrounded by additional 10 C 6 fragments. At the model potential level of theory, optimization returned the planar C 50 H 18 cluster. In contrast, PM3 returned two symmetrically equivalent curved global minima of C 2 symmetry, which correspond to global minima (E Rel curved = 0.0 kcal mol, no imaginary frequencies, with the first real frequency 31.97 cm −1 ) and planar transition state of C 2h symmetry with relative energy E Rel planar = 5.641 kcal mol, two imaginary frequences of 67.88i and 53.64i cm −1 , and the first real mode of 31.56 cm −1 . A small increase of the relative energy of planar C 50 H 18 cluster in comparison with C 40 H 16 one is caused by enhanced π-resonance due to increase of the number of surrounding π-conjugated C 6 rings from 8 (C 40 H 16 ) to 10 (C 50 H 18 ). The R-B3LYP/6-31G * relative energy of planar transition state of C 2h symmetry in respect to equivalent curved global minima conformers of C 2 symmetry is very small and equal to E Rel planar = 0.504 kcal mol. Comparison of the relative energies of planar C 2h transition states of C 40 H 16 and C 50 H 18 clusters clearly demonstrates increasing of π-resonance energy with increasing of the number of C 6 fragments in cluster backbone, which competes with deformation energy caused by introduction of 5/7 cores in the carbon lattices of graphene flakes. PM3 improper torsion angles along and perpendicular to 5/7 cores (figure SI1-1) for curved global C 50  An extended C 96 H 24 cluster was considered to study structural relaxation of phagraphene core with two 5/7 cores embedded into hexagon lattice environment (figure SI1-2). Two central 5/7 cores are separated by 3 C 6 rings and surrounded by additional external 30 C 6 fragments. Model potential optimization returns two equivalent global minima of C i or C s symmetry with relative energy E Rel curved = 0.0 kcal mol (no imaginary modes, the minimal frequency 37.23 cm −1 ), and planar C 2h transition state (relative energy E Rel flat = 2.10 kcal mol, one imaginary mode of 27.66i cm −1 , the first real frequency 31.37 cm −1 ) which is much smaller than the accuracy of the force field approach itself. Extension of C 50 -based cluster up to C 96 leads to prevalence of mechanical stress over π-resonance on the formation of the distinctive wave-like cluster even at the force field level of theory. The PM3 method returns a significantly larger potential barrier of E Rel flat = 44.33 kcal mol between the ground state minima (no imaginary modes, the first real mode is 17.33 cm −1 , E Rel curved = 0.0 kcal mol) and planar transition state, which has 3 imaginary modes of 126.78i, 77.22i, and 66.03i cm −1 with the first real one of 27.87i cm −1 .
Two planar C 32 H 14 clusters of C 2v symmetry (figure 4(c)) can be fused through the carbon atoms localized at C 7 external sides along 5-7 direction (figure 5(c)). One can expect that the fusion and completion of the lattice by additional 22 carbon may lead to formation of planar extended rhombus C 86 H 26 flake. The atomic structure of the flake was optimized by model force field, semiempirical PM3 and ab initio DFT B3LYP/6-31G * with both restricted and unrestricted M = 1 and 3 multiplicities. At all levels of theory, two symmetrically equivalent saddle shape conformers which correspond to equivalent global minima (figure 5(d)) of C 2v symmetry with relative energy E Rel curved = 0.0 kcal mol were revealed. The planar transition state of D 2h symmetry which separates the minima has MM+, PM3 and U-B3LYP/6-31G * relative energies E Rel planar = 15.85, 0.557 and 0.329 kcal mol/atom. The U-B3LYP/6-31G * open-shell triplet is the ground state of both planar and saddle conformers with singlet-triplet splitting equal to 13.708 kcal/mol for saddle conformer. The MM+ and PM3 Hessian calculations return groups of imaginary frequencies from 64.94i to 43.52i cm −1 and from 78.51i to 9.46i cm −1 , respectively, for planar transition state with the first real modes of 18.37 cm −1 and 13.76 cm −1 . Structural stress of the C 86 H 26 cluster is caused by the symmetrical mismatch of 5/7 cores with a perfect hexagonal lattice of graphene. Distortion of initial planar structure and high relative energies of planar transition state C 86 H 26 conformer at all 3 levels of theory in the range of 15.9-0.3 kcal/mol/atom unequivocally demonstrate strong non-linear response of the atomic structure for extension of the carbon lattice due to the local character of the π-resonance stabilization in comparison with completely delocalized mechanical stress over the whole flake.
Let us consider both planar and curved conformers of C 2v symmetry C 40 H 14 flake composed of 5/6/7 core (figure 6(a)) surrounded by 10 carbon hexagons. The PM3 planar conformer is open-shell singlet with relative energies of closed-shell singlet and triplet states equal to 24.078 and 2.196 kcal mol, respectively. The B3LYP/6-31G * revealed that the planar C 40 H 16 is a closed-shell singlet with relative energy of triplet state equal to 12.349 kcal/mol. At the PM3 level of theory, both singlets and triplet of the planar conformer have 4 imaginary vibrational frequencies (table SI1-1), making the conformer being a transition state between two symmetrically equivalent degenerated curved global minima (figure 6(a)).
The departure from C 40 H 16 planar C 2v transition state, to distinctive non-planar C s conformer (figure 6(a)) leads to visible energy gain of 37.254 kcal/mol for PM3 ground open-shell singlet and 17.408 kcal/mol for B3LYP/6-31G * ground closed-shell singlet. The relative energies of closed-shell singlet and triplet states at the PM3 level of theory are equal to 27.398 and 0.401 kcal mol, respectively, making open-shell singlet and triplet states almost degenerate with very small singlet-triplet splitting, which is smaller than the accuracy of the PM3 method. The B3LYP/6-31G * relative energy of the triplet state is equal to 6.667 kcal/mol. At the PM3 level of theory, no imaginary frequencies in the vibration spectrum of the non-planar C 40 H 16 cluster of C s symmetry for closed-shell singlet and triplet states are observed, making them the global and local minima at potential energy profile of C s C 40 H 16 conformer.
Two antiparallel 5/6/7 cores can be fused through 4 carbon hexagons with formation of extended sp 2 C 62 H 20 cluster (figure 6(b)), in which the double 5/6/7 core is surrounded by 12 hexagons. The PM3 open-shell singlet is energetically favourable than the closed-shell singlet and triplet states with relative energies of 40.034 and 3.164 kcal/mol, respectively. The PM3 vibration spectra of closed-and open-shell singlets revealed five and seven imaginary frequencies, respectively (table SI1-2), which indicate that the planar conformer is a transition state in both spin states. At the B3LYP/6-31G * level of theory, only a closed-shell singlet planar conformer was located.
Departure from sp 2 planar C 2h symmetry leads to localization (force field, PM3 and B3LYP/3-21G) of curved C 2 C 62 H 20 conformer of complex shape (figure 6(b)). Reduction of symmetry from PM3 planar C 2h transition state to two equivalent opened-shell singlet curved C 2 global PM3 minima leads to visible energy gain of 65.056 kcal/mol with relative energies of non-planar closed-shell singlet and open-shell triplet states of 49.245 and 13.157 kcal/mol, respectively. No negative frequencies in vibrational spectra for all 3 spin states of non-planar conformers were observed.
For all considered finite clusters, it was found that embedding of the fragments which contain pentagonand heptagon carbon rings in hexagonal environment in all cases leads to breakdown of perfect planar hexagonal symmetry of π-conjugated lattices with formation of curved flakes. Comparison of model potential and electronic structure calculations unequivocally demonstrates the leading role of mechanical stress in formation of the fine details of the lattices atomic structure. Following structural and symmetry analysis and using the Hessian calculations, it was found that planar conformers are just transition states with one or several imaginary vibration modes, which connect two symmetrically equivalent curved ground state minima on potential energy surfaces.  [53] which causes its transformation into dicyclopenta[ef , kl]heptalene (alternative names azupyrene or SW-pyrene) unit. It is necessary to note, that formation of SW fragments cannot be proceeded without distortion of both pentagonal and heptagonal azupyrene fragments and departure of the initial angle between translation vectors of graphene (P6/mmm space group) from 60 • to 74.02 • of phC(0,0).

Topological and quantum stability of 1D and 2D
The heptagon external angle of 127.53 • of a single freestanding azupyrene (figure 7(b)) is very close to one of freestanding regular planar heptagon (128.57 • ). The heptagon vertex of 7/5-5 junction is adjusted to two fused regular planar pentagons through the common zig-zag edge with the angle of 144.00 • . It is necessary to note, that for 2D phC(0,0) lattice this angle is equal to 141. It is necessary to note the dominant character of stretching force constant K r for sp 2 carbon atoms which is equal, for example, to 469 kcal/mol/Å 2 in respect to K θ bending sp 2 ∠C-C-C constant (85 kcal/mol/rad 2 ) for AMBER94 force field [100]. The dominant contribution of the stretching mode to the total energy leads to significant distortion of the heptagon angle (128.57 • for regular heptagon, see above) in 7/5-5 region of the planar 2D phC(0,0) lattice, making it 141.2 • , which is almost equal to the angle adjacent to 5-5 edge (141.3 • ).
For simplicity, the contribution of aromatic resonance to total energy can be ignored by effectively making the corresponding force constant responsible for keeping planar topology of any aromatic systems equal to 0. The accumulated structural mismatch (15.43 • -13.67 • ) of the units should lead to distortion of planar 2D phC(0,0) lattice with displacement of the atoms from the aromatic plane in the perpendicular direction to compensate the mechanical stress and to conserve the bond lengths minimizing the distortions of bending angles.
So, for low-dimensional cases (1D, 2D) a Topology Conservation Theorem (TCT) can be formulated in the following way: 'To conserve the planar topology of one-or few-atomic layer one-unit-cell-thick low-dimensional crystals with small or zero stabilizing force constants acting in the perpendicular direction to the lattice plain, and to avoid uncompensated mechanical stress, the free-standing constituting fragments (unit cells) must perfectly fit planar low-dimensional space. Due to the leading contribution of the stretching force constants to the total energy, mechanical stress accumulated by small regular structural mismatch of planar structural units should be compensated by out-of-plane structural deformations coupled with either multiplication of translation periods or breakdown of periodicity in one or two dimensions'.
Experimentally discovered perfect planar hexagonal and triangular lattices of graphene [9], h-BN [101], graphane [102], graphdiynes [103,104], gamma- [105] and holey-graphynes [106], g-C 3 N 4 [107], and g-C 4 N 3 [108] perfectly satisfy mandatory requirements of TCT since the structural fragments of low-dimensional crystals perfectly fill in the planar 2D space. It is necessary to note that hypothetic free-standing one-side hydrogenated graphan, so-called graphone [109], does not satisfy TCT mandatory requirements due to strong structural anisotropy. In fact, experimental graphone lattice is stabilized by Ni(111) surface [109] by weak van-der-Waals interactions and, probably, by rather energetic chemisorbtion forces induced by direct interactions of dangling carbon bonds with the substrate in interface region. Reversible dehydrogenation of graphone resulted in complex desorption processes with coverage-dependent changes in the activation energies for the associative desorption of hydrogen as molecular H 2 [109].
Two TCT corollaries can be formulated as well: first, any kind of external either positive or negative pressure may effectively stabilize planar 2D lattices prone for deformations by annihilation of the internal structural stress; second, in stochastic atomistic lattices which are formed by structural units which do not perfectly fill the low-dimensional space, the structural stress and consequent structural regular deformations cannot be accumulated because of mutual compensation of displacements in opposite directions.
The truthfulness of the first corollary can be illustrated by experimental fabrication of InGaAs/GaAs micro-and nanotubes with controllable inner diameter (from 4 μm to 4 nm) by selective etching of stabilizing sacrificial AlAs substrate [110]. The GaAs and InAs pair demonstrate 7.2% lattice structural mismatch, so after freeing of the (2ML)GaAs/(2ML)InAs vertical heterostructure by breakdown of stabilizing substrate, the heterostructure spontaneously rolls into a nanotube with GaAs fragment forming the inner layer of the nanotube. In this method, the AlAs substrate effectively provides stabilizing negative pressure which keeps the external GaAs/InAs vertical heterostructure perfectly planar during initial stages of the synthesys.
The truthfulness of the second corollary can be illustrated by the experimental synthesis of an sp 2 -hybridized one-atom-thick flat carbon membrane with random arrangement of four-, five, six and seven carbon polygons by electron irradiation of graphene [61,111]. Random distribution of the defects effectively eliminates periodicity of parent sp 2 graphene lattice and accumulation of mechanical stress caused by introduction of different types of regular defects.
The first region clearly indicates the extent of aromatic stabilization of π-system up to 12 unit cell length (up to ∼58 Å) with clear 1/x energy dependence. From 11 unit cell length (∼52 Å) the energy gain per unit cell demonstrates perfect linear dependence with an almost constant value close to −1 kcal/mol/unit cell. The 14 × 2 and 16 × 2 nanoribbons demonstrate some secondary aromatic antiresonances which slight decrease of the energy gain (∼0.03 kcal/mol/unit cell). The last 18 × 2 nanoribbon demonstrates a small positive deviation from perfect linear dependence due to the lack of accuracy of optimization. For the sake of clarity, the 14 × 2, 16 × 2 and 18 × 2 nanoribbons are presented as separate dots on the PM3 curve.
PM3 calculations of extended n × n (n = 3-8) phC(0,0) flakes (figure SI1-3) also demonstrate negative stabilization energy of screw conformers from −0.32 to −0.72 kcal/mol/unit cell without conclusive energy dependence. The absence of clear 1/x dependence (figure 7(f)) may be explained by the possible existence of numerous curved isomers among which the true global minimum is difficult to locate.
The phonon dispersion (figure 7(g)) of 2D planar SW lattice calculated at the PW DFT PBC level of theory reveals no imaginary modes. This contradicts to concerted previous results (structural analysis of phC(0,0) lattice, cluster calculations of finite n × 2 and n × n phC(0,0) nanoribbons and nanoflakes, Hessian calculations, which reveal imaginary modes for all planar conformers without imaginary modes for curved ones, analysis based on newly formulated TCT theorem) for finite phC(0,0) flakes presented above, which unambiguously demonstrate the transition state nature of planar 2D phC(0,0) conformers. This obvious failure of phonon dispersion calculations can be explained by three major reasons: first, TI is a symmetry restriction introduced by PBCs, which obviously leads to artificial stabilization of a perfect planar lattice. Second, to calculate the phonons with arbitrary (and large) k-vectors, proper potential should be used with a supercell 2-3 times larger than the real translation vector of the relaxed lattice. Since the true period of phC(0,0) is at least 32 unit cells (figures 7(d) and (e)), at least 64 × 64 supercell should be used to calculate phonon spectrum of the SW lattice. Third, and probably the most important reason is that the minimization of planar 2D lattice using the PBC approach was performed in complete active space of normal coordinates, which includes only subspace of normal coordinates which satisfy TI conditions. It is necessary to note that for periodic lattices which satisfy TI conditions, the dimensionality of TI SS NC coincides with CAS NC by definition. The mismatch of structural units of phC(0,0) leads to its departure from planar 2D lattice with formation of aperiodic (2D case) or superperiodic screw (1D case) conformers for which aperiodic coordinates should be included in CAS NC. It effectively makes CAS NC dimensionality much greater than TI SS NC one. From a mathematical point of view, minimization of phC(0,0) in TI SS NC is in fact constrained minimization which can artificially return a regular point on true CAS NC as a local or global minimum or maximum. Consequent phonon calculations using geometry obtained by constrained minimization procedure, as a result, may artificially return no imaginary modes in the dispersion law.
Summarizing the paragraph, it was found that structural units of phC(0,0) lattice do not fit each other and cannot perfectly fill in its planar low-dimensional crystalline space. Assuming the force constant perpendicular to phC(0,0) plane is negligibly small in combination with a dominant character of stretching force constants in comparison with bending ones, a TCT was formulated and proved. Two theorem corollaries were referenced and illustrated by known experimental results as well. It was shown that a mismatch of azupyrene structural units should lead to breakdown of phC(0,0) planar topology. In accordance with newly formulated TCT, the force field and electronic structure calculations of finite phC(0,0) nanoribbons and flakes resulted in bent and screwed structures as global minima with planar conformers between two symmetrically equivalent global minima as transition states. In particular, at the PM3 level of theory, it was found that the screw period of rotation of n × 2 phC(0,0) nanoribbon is equal to n = 32. At all levels of theory, Hessian calculations unequivocally confirmed transition state nature of all finite planar phC(0,0) nanoribbons and flakes revealing single or even multiple imaginary modes in their vibrational spectra. Phonon dispersion calculations of planar 2D phC(0,0) revealed no imaginary modes that was interpreted as a result of constrained minimization of phC(0,0) lattice in a subspace of normal coordinates which satisfy TI conditions. The constrained minimization in finite TI SS NC instead of Hilbert CAS NC leads to erroneous localization of a regular point in complete active space as a global/local minimum with the consequent calculation of phonon dispersion law in it. Since in this case the force matrix is calculated at a wrong point (not at a function extremum) the phonon dispersion law cannot be regarded as valid. It was concluded that for low-dimensional complex crystalline lattices with multiple sublattices phonon spectra calculations cannot be used as final and solid proof of structural stability in the case of broken periodicity. To avoid artificial stabilization caused by linear PBC symmetry restrictions, one must perform structural and symmetry analysis of a lattice first, considering symmetry and possible structural mismatches of structural units to satisfy TCT theorem mandatory requirements. Cluster calculations of finite fragments can be used to support TCT analysis as well.

phC(0,1) (phagraphene) lattice
As it was described above, to generate a regular phagraphene rectangular lattice [36], the constituting 5/7 fragments should be oriented in opposite directions, so structurally the neighboring 5/7 pairs cannot share the same vertices. The unrelaxed geometry of 5/7 pair is characterized by the 123.43 • external angle at 5/7 seam, which does not fit the 120 • of the neighboring hexagons. According to TCT, to avoid accumulation of internal structural stress, 2D and 1D phC(0,1) lattices should depart from perfectly planar topologies with the formation of a bent structures. To elucidate the symmetry effects and mechanisms of structural stabilization, planar and curved conformers of phC(0,1) finite narrow nA × mB nanoribbons (A and B are translation vectors along X and Y directions, respectively) with n = 1-25 and m = 1-18 unit cell length and width were studied using model potential MM+, semiempirical PM3, ab initio B3LYP/6-31G * , and plane wave PBE methods in cluster approximation. For the sake of comparison, the PBC PW PBE calculations of 1D and 2D phC(0,1) 1A × 1B and 1D 5A × 1B periodic supercells of wave (one wavelength) and planar phagraphene nanoribbons were used to reveal the effects of symmetry restrictions introduced by TI conditions of PBC. First, both one-unit cell 1D and 2D calculations resulted in perfectly planar lattices, whereas 5A × 1B periodic supercells revealed two conformers, namely (i) the ground state wave of 5 unit cell wavelength and (ii) the planar high-energy phagraphene nanoribbon.
Optimizations of nA × 1B nanoribbons for n = 1-25 at various levels of theory revealed two distinctly different types of conformers: the lowest in energy ground state wave conformer with the wavelength of 5 unit cells (or 33.844 Å at PM3 level of theory. Each single wave supercell contains 137 carbon atoms of total weight of 1644 Dalton) and planar conformers that are approximately from 20 (model potential) to 34 (PM3) kcal/mol/unit cell higher in energy of corresponding wave structures (figures 8(a) and (b)). Cluster PW PBE optimizations of 1A × 1B, 3A × 1B and 5A × 1B flakes returned curved (n = 1, 3) or wave conformers as the ground states as well. Both ab initio B3LYP/6-31G * and PBC PW PBE approaches reveal close relative energies of 23-25 kcal/mol/unit cell for planar structures as compared to wave-like ones. Using ab initio R-B3LYP/6-31G * and U-B3LYP/6-31G * calculations it was found that for phC(0,1) (n = 1-4) nanoribbons the triplet and closed-shell singlet states are almost degenerated in energy. For phC(0,1) 8A × 1B wave conformer, the open-shell singlet is the ground state with the energy of singlet-triplet splitting equal to 8.4 kcal/mol/unit cell. The energy gains per atom (1.3 × 10 0 -7.7 × 10 −1 kcal/mol/atom) caused by structure wave formation of phC(0,1) is comparable or greater than the energy of van der Waals interactions (9.6 × 10 −1 -9.6 × 10 −2 kcal/mol/atom, [112]). Perfect correspondence of model potential and semiempirical structural parameters and relative energies of all nA × 1B flakes with the results of cluster and PBC DFT calculations unequivocally indicates the structural origin of the waving effect. The smallest energy gain is revealed for the shortest 1A × 1B cluster which is equal from 2 (model potential) to 13 (PM3) kcal/mol/unit cell while B3LYP/6-31G * exhibits 6 kcal/mol (figure 8(c)). One can see the sharp drops in relative energy for half-(n = 3) and one-(n = 5) wavelength structures. Consistent elongation of the length of finite nanoribbons for n = 10 (maximum length for B3LYP/6-31G * approach), n = 20 (semiempirical PM3) and n = 25 unit cells (model potential) leads to concerted oscillations of the energy gain per unit cell with local minima correspondent to integer numbers of half-waves at n = 8, 10 (B3LYP/6-31G * , PM3, model potential), 13, (PM3, model potential), 18, 20, 23 and 25 (model potential). For PM3, the elongation of the nanoribbon beyond n = 14 leads to almost constant (33.4-33.6 kcal/mol/unit cell) energy differences between wave and planar conformers.
In table SI1-3 the model potential and PM3 vibration frequencies of finite nA × 1B phagraphene nanoribbons are presented. It was found that the accuracy of model potential structure optimization mostly determines the results of Hessian calculations and should be 10 −5 kcal/mol/Å or even better. The MM+ wave clusters except for n = 13, 20 and 25 have no imaginary frequencies in vibration spectra. Small imaginary frequencies (12i-5i cm −1 , table SI1-3) of MM+ n = 13, 20 and 25 finite nanoribbons are likely caused by either poor convergency during structure optimization or probably even some Hessian algorithm flaws. In contrast with the wave conformers, all n = 1-25 finite planar nanoribbons have several distinctive imaginary frequencies (76i-2i cm −1 , table SI1-3) in vibration spectra.
At PM3 level of theory, the n = 1-17 finite nA × 1B wave nanoribbons conformers have no imaginary frequencies in vibration spectra. Planar conformers with n = 1-18 have one or more imaginary frequencies (135i-71i, table SI1-3), which directly indicates that they are transition states between two equivalent ground state waves. It is necessary to note that the GAMESS code erroneously returns no imaginary modes for n = 12-15 and skips many of them for n = 10-18.
The phC(0,1) 25A × 2B and 15A × 2B nanoribbons calculated by both model potential and PM3 levels of theory keep the same, n = 5, wavelengths of structure waves. At the model potential level of theory, the nA × 4B (n = 2-72) phagraphene finite nanoribbons (figure 9, table 1) also reveal wave behavior with two localized waves of 18 (ground state, 0.0 kcal/mol/unit cell relative energy) and 6 (1.999 kcal/mol/unit cell) unit cell length, which can be considered as the second in energy excited state. Nanoribbons of 72 unit cell length correspond to 4 wavelengths of the ground wave cluster and 14 wavelengths of the short-wave conformer. The planar conformer with the highest relative energy of 2.618 kcal/mol/unit cell corresponds to the transition state between equivalent ground state waves of either 18 or 6 unit cell wavelengths.
At the model potential level of theory rectangular nA × nB (n = 2-18) phagraphene flakes (figure 10) also reveal the same wave behavior with the wavelength of 18 unit cells. For n = 8 almost planar structure is located, which is reflected in the half-wave structural antiresonance.
The dependence of the wavelength of phagraphene structure waves upon the width and length of finite nanoribbons, which is in fact the wave boundary conditions, localization of different conformers with different wavelengths and convergence of the ground state waves for different widths to the same wavelength  (18 unit cells), indicate that this phenomenon can be interpreted in terms of homogeneous D'Alembert wave equation: where Δ is Laplace operator, v is phase velocity, and solution u = {ϕ n } is a set of standing waves each of which satisfies preset boundary conditions. In terms of structure waves, the finite phagraphene flakes could exist since the solutions of the homogeneous D'Alembert wave equation can be found.
should be considered, where f(x, t) is a function of non-zero external acting force of one or two independent spatial coordinates. The solution of the inhomogeneous D'Alembert wave equation is an infinite superposition of plane waves each of which is a solution of D'Alembert equation. It is necessary to note that in the particular case of phC(0,1) lattice, the structural waves could be regarded as the analogs of acoustic phonons since the lattice atoms demonstrate coherent displacements from the plane without polarization effects (figures 8(a)-10). Infinite unconstrained averaged in time low-dimensional wave-shaped lattices correspond to completely undetermined lattice parameters and atomic coordinates because the plane-wave functions with different wave numbers and phases span the whole momentum and spatial spaces.
Each single structural wave of phC(0,1) nA × 1B is a solution of the wave equation with corresponding de Broglie wave length λ, eigenstate, and eigenvalue. The wave can be considered as a quantum quasi-particle with defined non-zero absolute value of the impulse |p| = h/λ (where h is the Planck constant) with two symmetrically equivalent opposite wave propagation directions which makes the average impulse p = 0. Disregarding the effect of spontaneous departure of free-standing graphene from perfectly planar lattice with formation of random intrinsic ripples in graphene perpendicular to the crystalline plane (see, for example, [115][116][117]), for pristine graphene lattice, which demonstrates the structural wave length λ = ∞, the absolute value of structural impulse |p| = h/∞ = 0, which makes 1D and 2D graphenes dynamically stable lattices. The quasiparticle estimation of structural wave energy of phC(0,1) lattice (see SI2 section) returnes negligibly small energy per carbon atom (6.97 × 10 −11 kcal · mol −1 · C atom −1 ), which makes the lattice dynamically stable.
For the infinite unconstrained averaged in time low-dimensional wave-shaped lattices the infinite superposition of the structural waves of different wavenumbers and phases directly leads to completely undetermined lattice parameters and atomic coordinates. Any kind of structural constraints and/or physical perturbations may lead to wave function of structure wave collapse with instantaneous localization of structural parameters and atomic positions of the atomistic lattice.
Disregarding thermal fluctuations, the energy of a perfect structural wave infinite periodic lattice is completely and exactly determined by its structural parameters like wave length and amplitude. Following the quantum uncertainty principle, the coordinate uncertainty Δx = /ΔE = ∞, where ΔE = 0 is Heisenberg energy uncertainty, which makes the atomic coordinates of the wave completely undetermined before any kind of measurement. During structural spectroscopic experiment the wave function of undisturbed free-standing infinite structural wave quasiparticle should collapse with localization of atomic positions at random coordinates with breakdown of lattice periodicity and formation of an aperiodic atomic lattice. Such type of structural instability can be regarded as a lattice quantum instability which can be defined as a lattice instability of wave-shaped low-dimensional infinite free-standing crystalline lattices caused by complete quantum uncertainty of atomic positions before any kind of structural experiment or another external impact which should lead to the collapse of the wave function of unconstrained and undisturbed structural wave quasi-particle at random atomic positions and lattice parameters.
Using quasiparticle approach (see SI2 section) and ideal gas laws, the internal pressures/temperatures of two localized infinite phC(0,1) structural waves of 5 unit cell wavelength (figure SI1-4) and 18 unit cell wavelength (figure SI1-5) were estimated to be 2466.7 atm/2870 K and 569.2 atm/2668 K, respectively. Based on the above estimations, one may conclude that the unconstrained and undisturbed infinite phC(0,1) lattice before the collapse of wavefunction of structural wave quasiparticle during structural experiment is characterized by completely undetermined structural parameters and random distribution of atomic positions with averaged carbon-carbon distances in the range of 2.5-4.0 Å, 2312.5-569.2 atm internal pressure and 2818-2668 K effective temperature. Such conditions seem to be enough [61,111] for fundamental structural transformations with breakdown of perfect wave structure and periodicity lift-off with the formation of random aperiodic carbon atomic lattices. It makes impossible the existence of unconstrained and undisturbed infinite phC(0,1) lattice due to lattice quantum instability. This conclusion can be extended for the case of screw-shaped 1D phC(0,0) lattice since regular screws demonstrate periodic behavior as well and can be considered as structure waves which prone to lattice quantum instability effects.
Planar 2D graphene lattice can be considered as a wave with infinite wavelength and 0 amplitude for which the volume of structure wave box is equal to 0. It makes the lattice parameters and atomic positions determined, the effective temperature and pressure of graphene lattice equal to 0 K and 0 atm, respectively, and averaged carbon-carbon distances equal to its solid-state limit of 1.5 Å.
The PW PBE PBC calculations of planar and wave conformers of narrow 1D phagraphene nanoribbon of one unit cell width were carried out for the supercell consisted of 5 unit cells ( figure 8(a)), which corresponds to one wavelength of the lowest in energy wave conformer. The crystallographic group of phagraphene is Pmg with the rectangular unit cell having a, b unit vectors equal to 8.095, 6.652 Å, respectively. It was found that the planar conformer is significantly higher in energy with relative energy equal to 24.762 kcal/mol/unit cell. For the sake of comparison, the calculation of the narrow 1D phagraphene nanoribbon of one unit cell width using just one unit cell as a slab was performed, which generated a planar 1D crystalline lattice as it was expected. The crystalline lattice of planar 2D phagraphene was calculated using the PW PBE PBC approach following the parameters of the original publication [36]. Previously [36,42] the phonon dispersion of phagraphene was calculated using both PHONON and GULP [46,47] codes for 2 × 2 and 3 × 3 [36,45] supercells and no imaginary frequencies were detected (figures SI1-6(a) and (b)). The phonon dispersions of planar 2D phagraphene were calculated for 1 × 1, 1 × 2, 2 × 1, 2 × 2, 2 × 3, 3 × 3, and 4 × 1 supercells (figures 11(a)-(f)). In contrast to previous publications [36,42] the phonon dispersion for all considered supercells is characterized by distinctive imaginary frequencies. Visualization of the imaginary frequencies demonstrates that they correspond to wave-like vibrational excitations of phagraphene crystalline lattice in the out-of-plane direction. The phonon dispersion of narrow planar 1D phagraphene 4 × 1 nanoribbon is characterized by imaginary frequencies as well (not presented in figure 11). One can speculate that obvious failure to determine imaginary frequencies by both codes [34,44,46] could be caused by 2 reasons. First, since the structure wavelength of 2D phagraphene (18 unit cells) is much larger than 2 × 2 and 3 × 3 supercells used for phonon calculations, the atomic vibration potentials used to calculate phonons were not appropriate to determine the imaginary modes. The second reason for the failure may be caused by possible π-resonance stabilization of extended 2D aromatic lattices with possible transformation of the planar transition state of 2D phagraphene into metastable intermediate conformer. In both cases, routine PBC DFT calculations of 2D phC(0,1) based on optimization of the structure using just one unit cell were failed to locate its true atomic structure determined by wave behavior of the lattice.
The results of ab initio NVT MD PW PBE PBC and cluster simulations of various 2D, 1D, and 0D phC(0,1) phagraphene lattices at 300 K are presented in figure 12. 5 ps AIMD NVT PBE PW PBC simulations of 2D phagraphene ( figure 12(a)) lead some departure of carbon atoms from the lattice plane with formation of small amplitude irregular wave-like structure. Due to NVT algorithm restrictions which froze the MD box dimensions to optimized PBC sepercell size, no one could expect the lattice can be transformed from the planar to a wave conformer since the formation of structural waves should lead to significant contraction of the lattice. For example, at ab initio DFT B3LYP/6-31G * level of theory, the length of typical phC(0,1) 10A × 1B flake contracts from 73.24 Å for planar conformer to 63.60 Å (almost 10 Å absolute difference) for the wave one, which corresponds to relative 13% contraction. During the whole course of 5 ps NVT PBC MD simulations, the total energy of 2D phagraphene lattice fluctuates (±1.25 eV) around thermostated equilibrium level of −2890 eV without any traces of phase transition.
The AIMD PBC simulations at 300 K of planar 1D phagraphene phC(0,1) nA × 1B nanoribbon ( figure 12(b)) return very close result with some thermal expansion of the lattice with small departures of the atoms from nanoribbon plane and formation of wave-like irregular conformer. Like in the case of AIMD PBC simulations of 2D phagraphene, the MD simulation box imposes strict limitations on the topology, volume of the nanoribbon slab. The length of the ribbon is equal to 73 Å, which is significantly longer the length of the wave conformer (64 Å) and, in fact, is equivalent to strong stretching force. The total energy of 1D phC(0,1) nA × 1B nanoribbon does not reveal any sign of phase transition fluctuating around ±1.0 eV of 300 K thermostated equilibrium level of −1269 eV.
The AIMD PBC simulations at 300 K of wave conformer of phC(0,1) nA × 1B nanoribbon (figure 12(d)) conserve its wave shape. Its thermostated equilibrium total energy approximately equal to −1273.5 eV which is significantly lower the thermostated equilibrium total energy of corresponding planar phC(0,1) nA × 1B nanoribbon conformer ( figure 12(b)) equal to −1269 eV. The fluctuations of the total energy of the wave conformer is much smaller and lie in the region of ±0.5 eV without any traces of phase transitions.
The AIMD simulations of 0D planar phagraphene phC(0,1) 10A × 1B C 248 H 56 flake (figure 12(e)) demonstrate completely different behavior of both total energy and the lattice structure of the flake upon simulation time. Average total energy of the flake at the beginning of simulation was equal to ∼−2448 eV. During the course of simulation the total energy plummeted up to −2455 eV (7 eV drop, ∼0.5 ps), rised to −2450 eV (5 eV rise, ∼1.2 ps), another drop to ∼−2456 eV (6 eV drop, ∼1.9 ps), and finally equilibrated at ∼−2453 eV (another 3 eV rise, ∼3 ps). Dramatic energy change during the course of simulations is coupled with planar topology breakdown and structural transformation of initial perfectly planar flake to irregular wave conformer, which can be interpreted in terms of typical phase transition with fundamental change of atomic lattice.
Overall, it was found that the phC(0,1) lattice forms structure waves with wavelengths varying from 5 to 18 unit cells depending on the boundary conditions with planar 1D and 2D conformers as transition states between two equivalent ground state waves. At least one second in energy structure wave of 6 unit cell wavelength was localized as well for nA × 4B nanoribbons. In contrast with previous publications [36,42], all Hessian and phonon dispersion calculations of phC(0,1) revealed imaginary phonon modes for planar conformers. It was shown that NVT AIMD PBC algorithm artificially confines periodic lattices inside MD simulation box of constant dimensions which prevents their thermal relaxation and equilibration. In fact, NVT AIMD PBC approach effectively imposes artificial stretching or compressing forces keeping the dimensions of periodic slabs constant during the whole course of MD simulations. It was found that infinite 1D and 2D phC(0,1) lattices suffer quantum instability, whereas finite cluster or confined lattice of phC(0,1) may exist and could be experimentally synthesized.
In addition to 5/7/6 fused fragments, the structural mismatch for which is equal to 3.43 • , the phC(1,0) lattice is characterized by 5/5/7 and 7/7/5 fragments (figures SI1-7 and 13(a)) as well. The unrelaxed geometry of 5/5 pair is characterized by 144 • external angle at 5/5 seam, which is far beyond the 128.57 • of the perfect heptagon. The 7/7 pair has 102.86 • which does not match 108 • of the perfect pentagon. So, following TCT, the phC(1,0) should suffer internal structural stress, which should lead the lattice to depart from the perfect plane with formation of a bent structure.
The atomic structure and mechanical stress of finite phC(1,0) nanoribbons, flakes and nanotubes were studied by model potential, semiempirical PM3, ab initio B3LYP/3-21G and PW PBE methods. For the sake of comparison, PBC PW PBE calculations of 2D phC(1,0) lattice were performed as well. The phC(1,0) bending energies per unit cell of nA × 1B phC(1,0) nanoribbons (figure 13(a)) calculated as where E b is the binding energy of finite nA × 1B phC(1,0) nanoribbon, E cand E p are the total energies of curved and planar conformers, respectively, and n is the length of the nanoribbons in unit cells. At all levels of theory optimizations of one unit cell 1A × 1B phC(1,0) fragment returned planar atomic lattice (top-left of figure 13(a)).
The model potential optimization of finite fragments of nA × 1B phC(1,0) nanoribbons (n = 1 − 14) returned curved atomic lattices for (n = 3 − 14) which may be finally rolled to a multiwalled nanotube (an example of 12A × 1B phC(1,0) fragment is presented in figure 13(a)) or fused into a short single-wall phC(1,0) (14,0) zig-zag nanotube of 1 unit cell length. At the model potential level of theory, the bending energy varies from 0 kcal/mol/unit cell for n = 1, 2 to −12.53 kcal/mol/unit cell for n = 13 with distinctive 1/x behavior for (n = 3 − 11). A relatively sharp drop of E b for (n = 12, 13) is caused by the attractive interaction of both ends of the nanoribbon caused by its rolling. Following its universal curvature, the edges of phC(1,0) nanoribbon can be fused to form the (14,0) zig-zag nanotube of 21.368 Å diameter. The nanotube bending energy (blue dot in figure 13(a)) belongs to the same 1/x (n = 3 − 11) pattern due to the absence of vdW interactions between free edges and rather large nanotube diameter (and the number of constituting unit cells) to eliminate structural differences between complete and uncomplete nanotubes.
In contrast to curved conformers, all planar finite-length nanoribbons, except n = 1, 2 ones for model potential and n = 1 for PM3 and ab initio B3LYP/3-21G, demonstrate multiple (for large n even more than 50) imaginary modes, which makes planar conformers the transition states between two symmetrically equivalent minima.
It is necessary to note two very important features of vibration spectra calculations for one-atom-thick extended nanoclusters: (i) to make solid conclusions about stability of one-atom-thick nanoclusters, the quality of optimization must be very high with the gradient <1 × 10 −5 kcal/mol/Å; (ii) for one atom thick extended clusters with the number of atoms more than 54, existing algorithms of calculations of vibration spectra are not stable. To avoid artificial absence of imaginary modes in vibration spectra of the clusters, one has to use at least 2 independent codes for calculations of molecular vibrations (in this work GAMESS and HyperChem codes were used to cross control the results), combined with TCT structural analysis.
The PM3 vibration spectra calculations of phC(1,0) (13,1) chiral nanotubes constituted by 16-23 unit cells ( figure 13(b)) revealed no imaginary frequencies, which unequivocally demonstrate their structural stability. The shortest nanotubes of 14 and 15 unit cells length revealed one small imaginary mode each of 4.68i and 2.17i cm −1 , respectively, which is caused by soft vibrational modes of the nanotube fragments in the vicinity of fusion regions, for which even enhanced accuracy of optimization maybe insufficient to calculate vibration spectra.
One can speculate that phC(1,0) can spontaneously form roll nanotubes rather than planar 2D lattice or (13,1) seamless chiral nanotube. For the sake of comparison the phC(1,0) (13,1) seamless chiral nanotube is presented in figure 13(b) as well. The phC(1,0) 15A × 4B roll was optimized at the model potential level of theory and it is clearly seen that weak vdW interactions form roll structure with curvature even greater than the nanotube.
PW PBE PBC calculations of finite 4A × 4B planar and curved phC(1,0) nanoflakes returned −0.288 kcal/mol/C atom energy gain for curved conformer, which is very close to PM3 result of −0.238 kcal/mol/C atom for 4A × 4B flakes (see above). PW PBE PBC calculations of phonon frequencies of 2D phC(1,0) crystalline lattice returned 1 imaginary mode ( figure 13(c)) which corresponds to a bending mode perpendicular to the lattice plane that is in perfect agreement with vibration spectra calculations of all planar finite fragments. Therefore, one can conclude that planar 2D phC(1,0) is structurally unstable.
Overall, the TCT analysis and computer simulations of 1D and 2D phC(1,0) lattices revealed that: (i) structural mismatch of phC(1,0) units leads to mechanical stress with consequent formation of 1D helices, (13,1) chiral, or rolled nanotubes; (ii) at all levels of theory curved structures are symmetrically equivalent global minima on potential energy surface separated by planar conformers as transition states; (iii) stabilization bending energy of phC(1,0) lattice is comparable or greater the vdW interactions with any kind of support; (iv) calculations of vibrational spectra either for finite clusters or 2D lattices require incredibly high (up to 1 × 10 −5 kcal/mol/Å gradient) convergency criteria for global minimum optimization. Curved global minima have no imaginary vibrational frequencies at all levels of theory.
(v) Phonon spectra calculation of 2D phC(1,0) lattice returns one imaginary mode; (vi) calculations of vibrational spectra of finite planar clusters are prone to obvious mistakes which contradicts either each other or TCT analysis.

phC(1,1) lattice
In figures 14(a) and SI1-10 a 4A × 4B fragment of a planar phC(1,1) flake is presented. Each single unit cell of phC(1,1) consists of two antiparallel 5/6/7 units separated by 4 fused C 6 fragments arranged in stripes and separated from each other by infinite lines of fused carbon hexagons.
The lattice is formed by fused 5/6 and 7/6 fragments both of which confine hexagons at the seam angles. Like in the cases of phC(0,0), phC(0,1), and phC(0,1), the phC(1,1) lattice does not satisfy the TCT requirements. In particular, the angles formed by fusion of pentagons and heptagons with hexagons are equal to 132 • and 111.43 • , respectively, which is far beyond 120 • of perfect hexagons. The phC(1,1) lattice cannot be perfectly filled in by undeformed pentagonal, hexagonal and heptagonal fragments and should be relaxed with the formation of a bent structure even its structural parameters (for example, the model potential a 1 = 12.935Å, a 2 = 8.013Å, α = 119.317 • and β = 60.683 • ) are rather close to a perfect hexagonal lattice.
Atomic structure, energetic characteristics and vibrational frequencies of finite planar and wave nA × mB n = 0.5 − x, m = 1 − y phC(1,1) finite flakes and nanoribbons were calculated using model potential, PM3 and B3LYP/3-21G methods. It was found that at all levels of theory planar 2D conformers (figures 14(a) and S10) are transition states between two equivalent wave structures (figures 14(b)-(d)) with structure waves propagating along fused hexagonal lines. In particular, the model potential 11A × 8B phC(1,1) flake (C 3270 H 150 , figure 14(c)) is characterized by bending stabilization energy equal to −32.574 kcal/mol/unit cell or −0.877 kcal/mol/C atom. The wavelength of structure nA × 1B phC(1,1) wave ( figure 14(b)) is equal to 4 unit cells that is even shorter than the wavelength of phC(0,1) nA × 1B lattices (5 unit cells) and does not depend upon the width of the ribbon, which can be attributed to confinement effect of lines of hexagons.
All PM3-calculated planar conformers of nA × 1B (n = 1-4) revealed 2, 4, 7 and 11 imaginary vibrational frequencies, which directly indicate their transition state nature. No imaginary modes were detected for any wave conformers. B3LYP/6-31G * minimization of the smallest 1A × 1B cluster returned only planar conformer without imaginary modes in vibration spectrum (for larger clusters nA × 1B (n = 2, 3) vibration spectra were not calculated). At PM3 and B3LYP/6-31G * levels of theory stabilization energy of formation of 3A × 1B wave conformer in respect to the planar transition state is equal to −1.447 and −1.170 kcal/mol/C atom, respectively, which is close enough to −0.877 kcal/mol/C atom model potential stabilization energy for 11A × 8B phC(1,1) flake. The 4A × 1B cluster ( figure 14(b)) was studied only at the PM3 level of theory, for which stabilization bending energy is equal to −1.716 kcal/mol/C atom.
Overall, it was found that the planar 2D phC(1,1) lattice does not satisfy TCT requirements which makes it topologically unstable suffering structural deformations perpendicular to the main lattice plane. Using model potential, PM3 semiempirical, and ab initio B3LYP/3-21G * simulations it was found that phC(1,1) lattice forms structure waves which propagate along narrow zig-zag ribbons formed by fused carbon hexagons. The ribbons confine the wave behavior of the lattice keeping the wavelength of phC(1,1) structure waves constant and equal to 4 unit cells. Stabilization bending energy of phC (1,1) lattice is estimated to be 1-2 kcal/mol/C atom, which is greater than the vdW energy, so planar 2D phC(1,1) lattice cannot be stabilized by any kind of supports. Because of structure wave behavior, an infinite free-standing 2D phC (1,1) sheet should be prone to quantum instability effects.

Penta-graphene topological instability
The topological and structural instability of penta-graphene [35] was recently studied [37][38][39] and it was shown that at the one-electron level of theory its planar conformer is just a regular point on the potential energy curve [37] between two equivalent non-periodic saddle structures. The penta-graphene lattice instability is entirely caused by structural reasons and it should be thoroughly investigated to reveal the true nature of its atomic lattice.
The sp 3 central sublattice of penta-graphene is characterized by strong departure of bonding angles (99.262-100.800 • ) from a perfect undisturbed angle of 109.471 • formed between sp 3 atoms. Top and bottom carbon dimer sublattices of penta-graphene also demonstrate strong departure of bonding (110.974-112.601 • ) and torsion (127.168 • ) angles from perfect sp 2 environment (120 • and 180 • , respectively), which means the penta-graphene 2D planar lattice does not satisfy mandatory TCT requirements and should suffer strong structural distortion due to internal mechanical stress.
Both planar and curved conformers of C 72 H 28 (3 × 3) and C 120 H 36 (4 × 4) (figure 15(a)) penta-graphene clusters were calculated at PM3 and ab initio B3LYP/6-31G * levels of theory. C 72 H 28 and C 120 H 36 planar conformers were minimized using symmetry restrictions keeping constant z-coordinates of central sp 3 carbon sublattice and relaxing x and y coordinates of all C atoms at particular z-coordinates of non-equivalent top and bottom sp 2 carbon sublattices. Global minima for planar clusters were minimized at z = ±0.60555 and z = ±0.60731 Å for carbon atoms which belong to top and bottom dimer sublattices for (3 × 3) C 72 H 28 and (4 × 4) C 120 H 36 clusters, respectively. It is necessary to note, that the thickness of both clusters (1.211 and 1.215 Å, respectively) just coincide with high accuracy with PBC results [35].
To evaluate the curvature of saddle conformers of penta-graphene narrow nanoribbons cut out in a and b directions (figures SI1-11(g) and (l), respectively) of two unit cell width were calculated using PM3 method. Both type of nanoribbons spontaneously form short fragments of nanotubes (which can be also considered as rings) of either 5.98/7.14 Å (the ring circumference of 8 unit cells), or 9.94/11.19 Å (the ring circumference of 17 unit cells) of inner/outer radii r i . It is necessary to note that the main axis of a nanotube formed by a-oriented nanoribbon is parallel to b-dimension and the main axis of a nanotube formed by b-oriented nanoribbon is parallel to a-direction, forming b-, and a-nanotubes, respectively. Corresponding curvatures κ i = 1 r i and the total curvature of the saddle point K = κ bi κ ao = −0.015Å −2 , where b i and a o indexes correspond to the inner radius of b-oriented ring and outer radius of a-oriented ring, respectively. To prove the stability of the rings, unlocked b and a rings (figures SI1-11(f) and (k), respectively) were calculated as well and it was shown that both rings keep their perfect round shapes and radius values.
The ratios of inner and outer radii of b-and a-rings are pretty close to each other and are equal to 0.602 and 0.638, respectively. Different radii of b-and a-oriented nanotubes indicate significantly different bending force constants associated with either b-, or a-directions of one-unit cell thick penta-graphene, making the lattice significantly anisotropic. Following the saddle shape of the penta-graphene flakes (figures 15(a), SI1-11) and taking into account significantly different curvature radii of b-and a-nanotubes, one can conclude that b-and a-associated bending force constants have opposite directions and different absolute values. It is necessary to note that since the inner and outer radii of a-oriented ring are significantly larger the b one, corresponding bending force constant should be visibly smaller.
Elongation of b-oriented ring (figure SI1-11(g)) up to 5 unit cells generates perfect nanotube with the same inner/outer radii r i of 5.98/7.14 Å. Unlocking the b-nanotube (figure SI1-11(h)) resulted in some increase of the effective radius of the tube due to competition and partial mutual compensation of the b-and a-bending force constants.
One can speculate that during a hypothetical growth of one unit cell-thick penta graphene lattice, fusion of opposite lattice edges is hardly possible due to kinetic reasons and misorientation of opposite edges. Instead of the nanotubes and rings, formation of corresponding b-and a-oriented 1D rolls (figures SI1-11(j) and (m)) is very plausible. Since the bending energy of one-unit cell thin penta-graphene lattice is approximately as twice as large as vdW binding energy, one could speculate that after achievement of certain radius at which the resulting curvature of external penta-graphene roll would be rather close to planar structure, strong bending stress should form secondary roll structure with structural passivation of the growing edge, preventing formation of large-diameter penta-graphene rolls.
Considering sp 2 one-atom thick lattices, the negatively curved surfaces may form several regular 0D, 1D or 3D solids of regular structure. 0D species contain quasi-spherical round giant and star-shaped closed-shell fullerenes and tori (see, for example, [118,119]), which have both negatively and positively curved regions. 1D regular species [118], namely corkscrew-shaped nanotubes also demonstrate both negatively and positively curved fragments. From mathematical point of view, the penta-graphene negatively curved lattice cannot form neither 0D or 1D closed-shell solids due to two reasons: first, small penta-graphene flakes (figure 15(a)) do not have positively-curved fragments; second, both 0D tori and 1D corkscrew-shaped nanotubes have non-zero thick minor radi, so the length of extra torus/corkscrew fragments of negative curvature should be longer the inner positively curved ones which is not possible without introduction/cutout of extra three-atom thick unit cells from external/internal positively/negatively curved fragments of the lattices.
3D triply periodic minimal surfaces (TPMS) were mathematically introduced and studied by Schwarz [85] and Schoen [86,120] long before discovery of any complex forms of nanocarbons back in 1890 and 1970. A number of regular 3D sp 2 -carbon TPMS schwarzites or buckygyms consisted of just negatively curved lattices were introduced (see, for example, [88][89][90]). Considering the TPMSes and carbon foams on the base of them, Hyde et al [121] mentioned that 'we find we have entered an enormously complex, and virtually unexplored, area of geometry'.
All schwarzites which were developed on the base of TPMSs have 3D periodic structures (see, for examples, figure SI1-12) with 0 thickness of negatively curved continuous non-crossing surfaces constituted by hexagons and larger polygons (heptagons, octagons and nonagons) without any residual mechanical stress. The last property implies that the force constants normal to any TPMS surface can be ignored. Zero thickness of the surface also means that the surface integrals of both sides of the surface should be symmetrically equivalent. The penta-graphene lattice do not mathematically satisfy both TPMS properties since it has non-zero thickness with strongly anisotropic force constant field.
Extended penta-graphene na × 8b (n = 5-8) flakes (figures SI1-11(a)-(e)) demonstrate complex structural behavior upon the flake composition. Initially, the 5a × 8b flake ( figure SI1-11(b)) demonstrate distinguished a-related curvature. An extension of the lattice in a direction (figures SI1-11(c)-(e)) leads to graduate mutual compensations of the force constants ( figure SI1-11(a)) oriented in opposite direction with cancellation of overall negative curvature of the na × mb flakes with formation of b-oriented linear 1D penta-graphene nanotube (figure SI1-11(i)) considered above. It is necessary to note that since both types of relaxed minimal energy penta-graphene nanotubes are highly symmetrical, they cannot demonstrate any kind of curvature in any direction.
Model In fact, formation of the secondary coil should prevent the lattice growth due to passivation of the nanoribbon edge which in turn should limit the maximum diameter of penta-graphene rolls from 26 to 65 Å. Since the secondary coil introduces strong lattice anisotropy, it should generate non-zero bending force perpendicular the main b-oriented roll axis, which should lead to complete breakdown of the lattice periodicity in any direction.
Phonon spectrum of planar 2D penta-graphene (figure 15(c)) calculated using PBC PBE PAW approach reveals a small (less than 1i cm −1 ) imaginary mode around Γ-point, that is in contrast to [35] where no imaginary frequencies were revealed. The imaginary part has a gentle slopping and is associated with out-of-plane vibrations that are normally presented in various 2D materials. Therefore such negative part around Γ-point may be considered as numerical error and ignored with considering a low-dimensional lattice as stable according to the PBC approach.
Infinite multidimensional complete active space of normal coordinates of penta-graphene includes both finite TI and infinite non-TI subspaces. From a mathematical point of view, PBC constrains energy minimization of the penta-graphene lattice to finite TI subspace of normal coordinates, which leads to artificial localization of a regular point in true CAS NC as the structural global minimum. Consequent phonon dispersion calculations using the geometry obtained by constrained minimization procedure may artificially return no imaginary modes witnessing false stability of the lattice.
Within symmetrical TI constraints, 2D planar penta-graphene [35] consists of 3 non-equivalent sublattices with one central sp 3 atom and two (top and bottom) sp 2 carbon dimer sublattices oriented perpendicular to each other. Since non-equivalent sp 2 sublattices are separated from each other by 1.2 Å, they generate uncompensated torque (see section 3), which as a consequence leads to dramatic structural distortion of the lattice. In fact, the linear PBC approximation is a symmetry restriction that forbids (if applied) folding of any kind of low-dimensional lattices especially for non-periodic bend structures or structure waves with a wavelength greater than one unit cell or, in general, applied supercells. One should take into account that the linear PBC approach leads to artificial stabilization of low-dimensional crystalline lattices which may suffer a complete breakdown of periodicity due to, for example, formation of curved structures, with consequent artificial luck of imaginary modes in phonon spectrum.
Overall, it was shown that: (i) mathematically one unit cell thin penta-graphene lattice can not exist in either perfectly planar 2D lattice, or 0D closed-shell cages like spherical fullerenes, carbon tori or 1D corkscrew coils; (ii) penta-graphene lattice can form either saddle-shaped 0D flakes of negative curvature of K = 0.015 Å −2 , or 1D perfect linear minimal energy 1D nanotubes of either 5.98/7.14 Å or 9.94/11.19 Å inner/outer diameters or rolls of finite diameters and complex nature; (iii) due to strong bending force constants perpendicular to the lattice, the rolls have very asymmetric structure due to formation of secondary coil which should passivate the lattice growth, limits external roll diameter up to 65-70Å, introduce strong lattice anisotropy and generate non-zero bending along the main axis, which should, in turn, completely eliminate periodicity of the rolls in any direction; (iv) PBC approximation causes artificial topological stabilization of 2D penta-graphene due to restrictions imposed by linear translation symmetry; (v) for small flakes, two perpendicular sp 2 carbon dimer sublattices of penta-graphene separated by 1.2 Å cause uncompensated torque which breaks down topology and periodicity with the formation of aperiodic saddle structures; (vi) the stress energy of carbon atoms in planar flakes with preserved coordination is twice as great the vdW energy, making impossible to stabilize planar penta-graphene clusters at any kind of supports; (vii) constrained TI SS NC energy minimization of penta-graphene artificially returns a regular point in CAS NC surface as the global minimum; (viii) phonon dispersion law calculations of penta-graphene at constrained minimum artificially return no imaginary modes, which compromise the phonon spectrum calculations as a solid and final proof of stability of low-dimensional crystalline lattices; (ix) small bent n × m (n, m 5) penta-graphene clusters or nanorolls could exist and in principle may be synthesized.

Topological instability of ladderane lattice
Originally, the ladderane lattice [65] was designed following cyclobutane motif ( figure 16(a)). The Pmma orthorhombic unit cell of ladderane consists of two non-equivalent carbon atoms. The 2D crystalline structure was calculated using the PW PBE PBC approach and it was found that C-C bond lengths are equal to 1.640 and 1.511 Å, respectively, with α = β = γ = 90, which is rather close to the original data of [65] (1.650 and 1.513 Å with α = β = γ = 90). Four-fold coordination of each carbon atom neither belongs to sp 3 -nor sp 2 -hybridization, so the lattice cannot be perfectly fitted either by tetrahedral sp 3 or triangular sp 2 fragments due to symmetrical reasons. One can conclude, that 2D planar ladderane does not satisfy the mandatory TCP requirements and could be artificially stabilized by the PBC approach. Phonon dispersion of ladderane ( figure 16(b)) reveals no imaginary modes in the whole momentum space, which perfectly coincides with original results [65].
To check topological stability of the lattice, minimization of two finite unsaturated C 54 and saturated C 54 H 30 ladderane 4 × 5 clusters was performed using the PW PBE approach (figure 16(c)). Minimization of unsaturated C 54 cluster revealed dramatic restructuration of ladderane lattice converting initial 4 × 5 buckled ladderane C 54 cluster of C 2v symmetry to perfectly planar D 2h hexagonal lattice (figure 16(c) left). Complete restructuration followed by symmetry breakdown of the initial 2D ladderane lattice may lead either to formation of 2D graphene or separate 0D graphene flakes, that may lift the long-range crystalline order of regular covalent bonds between the atoms spanned all over the crystal. Like in the case of penta-graphene, PBC artificially stabilizes rectangular buckled ladderane lattice keeping its shape and translational symmetry invariance. PW PBE energy minimization of saturated C 54 H 30 cluster (figure 16(c) right) leads to a complete breakdown of its symmetry, long-and even short-range orders with the formation of an irregular exotic atomic lattice with 3-, 4-, 5-, 6-, 7-, and even 8-member rings. Relaxation of both C 54 and C 54 H 30 clusters leads to complete decomposition of finite regular 2D rectangular ladderane lattice which satisfies TI conditions with the formation of either graphene flakes or irregular 0D finite carbon clusters with breakdown of the topology and both long-and short-range orders.
The internal forces of ladderane are strong enough to overcome the C-C bonds regrouping potential barriers and completely disintegrate the lattice with breakdown of translation symmetry and topology of the regular crystalline lattice. By definition, the phonon dispersion calculations are performed for the lattices which satisfy the TI conditions of TI subspace of normal coordinates x which belongs to Hilbert complete active space of normal coordinates, x ∈ X. In the case of ladderane, the X space of infinite dimension includes aperiodic normal coordinates of disintegrated fragments of ladderane as well as finite TI x, X x, of regular ladderane lattice. Because of it, the constrained extremum localized in TI x subspace is just a regular point on Hilbert potential energy surface located far away from any true global minimum/ maximum. Consequent calculations of phonon dispersion law at constrained TI SS NC minimum could artificially return no imaginary modes witnessing false stability of the lattice.
Summarizing the paragraph, it was found that 2D ladderane lattice does not satisfy TCT mandatory restrictions due to dramatic departure of the elementary ladderane rectangular fragments either from triangular (sp 2 carbon), or tetrahedral (sp 3 carbon) structural units. Linear PBCs lead to artificial stabilization of ladderane keeping intact its rectangular lattice with perfect translation symmetry. The lift of translation symmetry restrictions leads to complete breakdown of the integrity, translational and local symmetries of the lattice with formation of localized fragments caused by strong internal structural stress. In PBCs, the energy minimization of the ladderane 2D lattice is reduced to constrained minimization in subspace of translationally invariant normal coordinates far away from any true global/local minimum/maximum of Hilbert complete space of normal coordinates, which leads to artificial localization of a regular point on true potential energy surface as a special point for which no imaginary phonon modes were detected. 2D ladderane lattice suffers severe topological instability which cannot be detected by regular phonon calculations due to strong departure of 2D rectangular lattice from a true global minimum or maximum.

The mechanism of structural stabilization of 2D biphenylene lattice
All previous calculations of biphenylene lattice were performed for perfectly planar lattices applying symmetry restrictions either in cluster approach by D 2h symmetry point group [27], or PBC [28,122] Pmmm space group with six basis carbon atoms in the rectangular unit cell. The lattice constants of 2D biphenylene are a = 3.76 Å and b = 4.52 Å with 3 different C-C bonds, with the lengths of l 1 , l 2 , and l 3 equal to 1.45, 1.46, and 1.41 Å, respectively. No imaginary phonon modes were detected in previous phonon dispersion calculations [122].
To analyze topological stability of 2D biphenylene lattice (figure 17(a)) its parent constituting structural units, namely cyclobutadiene C 4 H 4 (figure 17(b)) cyclooctatetraene C 8 H 8 (figures 17(c) and (d)) and biphenylene (C 6 H 4 ) 2 (figure 17(e)) should be considered first. Molecular structure of cyclobutadiene, cyclooctatetraene both ground state tub and planar transition state conformers and biphenylene were calculated using ab initio B3LYP/6-31G * and PW PBE approaches. Cyclobutadiene C 4 H 4 belongs to [n]-annulene family (actually, it is the smallest [4]-annulene) and is characterized by its antiaromatic π-system with four π-electrons. It is a correlated molecule in which the Jahn-Teller effect distorts the carbon 4-member carbon cycle converting triplet carbon square to singlet carbon rectangule with two non-equivalent C-C bonds [123]. At the MR-CCSD(T) level of theory the ground state rectangular singlet is characterized by 1.483 and 1.291 Å long and short C-C bond lengths, respectively, whereas exited triplet square conformer has four equivalent 1.386 Å C-C bonds [124].
Antiaromatic cyclooctatetraene C 8 H 8 ( figure 17(d)) is another member of [n]-annulene ([8-annulene]) family. Its ground state tub-shaped conformer [125,126] is characterized by two C-C-C non-equivalent angles of 121.1 • and 117.6 • . Anionic form [C 8 H 8 ] −2 is an aromatic molecule with a perfectly planar structure of D 8h symmetry with eight ∠C-C-C equal to 135 • . The neutral transition state conformer C 8 H 8 , cyclooctatetraene-TS ( figure 17(c)), has a planar structure as well with reduced symmetry from D 8h point group to D 4h caused by the Jahn-Teller effect.
The free-standing C 12 H 8 biphenylene unit ( figure 17(e)) is a polycyclic molecule [127,128] with two aromatic benzene rings connected through strongly antiaromatic central rectangular cyclobutadiene fragment [129]. Since the central antiaromatic C 4 unit keeps its planar structure due to structural limitations, the whole molecule also keeps its planar structure as well. In the unit cell of periodic lattice, every single C 4 fragment of 2D biphenylene displays an almost square structure with two types of C-C bonds, equal to 1.4572 and 1.4536 Å, respectively, which is very close to PW PBE bond length of 1.444 Å of C 4 H 4 cyclobutadiene.
Every single C 6 hexagonal fragment in 2D biphenylene lattice has two types of ∠C-C-C angles equal to 110.0 • (two angles) and 125.0 • (four angles), which do not match neither perfect aromatic C 6 H 6 (120.0 • ) In 2D biphenylene, the most dramatic changes are revealed in C 8 fragments, which are very different from either cyclooctatetraene ground state non-planar D 2d tub conformer, or transition state planar D 4h one. C 8 of 2D biphenylene has 3 types of C-C bonds of 1.445 Å (2 bonds), 1.454 Å (2 bonds), and 1.405 Å (2 bonds) with 2 types of ∠C-C-C angles equal to 125.0 • (4 angles) and 145.0 • (4 angles), respectively, whereas tub conformer has non-zero torsion angle of 54.7 • , 8 equivalent ∠C-C-C angles of 127.3 • and two types of C-C bonds of 1.346 (4 bonds) and 1.469 Å (4 bonds). Planar TS conformer has 8 equivalent ∠C-C-C angles of 135.0 • and two types of C-C bonds of 1.350 (4 bonds) and 1.472 Å (4 bonds).
Overall, that 2D biphenylene lattice does not match all cyclobutadiene C 4 , benzene C 6 , cyclooctatetraene C 8 , or biphenylene C 12 H 8 conformers. It could be considered as a correlated lattice for which different types of correlations, aromatic and antiaromatic resonances may cause dramatic structural distortions of constituting fragments, the lattice symmetry and topology.
The structure and topology of different 2D biphenylene fragments like n × n flakes and n × m ribbons were studied using model potential, PM3 semiempirical (figure SI1-13) and ab initio B3LYP/6-31G * methods. Stabilization bending energies of the fragments calculated by PM3 semiempirical and two types of model potentials with (MM+(a)) and without (MM+(b)) taking into account the aromatic resonance [66,130]) are presented in ( figure 17(f)).
All dependencies of PM3 bending stabilization energies upon dimension of n × 5 and n × 4 nanoribbons and n × n flakes of biphenylene lattice (figure 17(f)) demonstrate clear 1/x patterns. It is necessary to note that some deviations from 1/x pattern for n 14 for n × 5 nanoribbon is caused by a lack of accuracy of structural minimization of extended ribbons with the number of atoms in the lattices >450. The energy gain of square-shaped flakes converged to −0.006 kcal/mol/C atom which is smaller than the experimental values of vdW energy [112]. PM3 structural optimization of n × 4 and n × 5 nanoribbons reveals −0.01 and −0.02 kcal/mol/C atom stabilization energies, respectively, which is comparable with the lowest limit of vdW binding. Structural minimization of n × 4 and n × 5 nanoribbons revealed structural curvature with radii of R n×3 = 87.266 Å and R n×4 = 85.139 Å, respectively ( figure SI1-13). The curvature radii of 42.761 × 36.287 × 4.042 Å 10 × 10 flake in x, y and diagonal x + y directions are found to be R x n×n = 172.5Å, R y n×n = 200.5Å, and R xy n×n = 167.8Å, respectively, with a maximum deviation from the plane equal to 4.042Å.
The MM+ algorithm [130] allows one to separate different stretching and torsion components to the total energy of carbon-based atomic lattices with active π-system. Let us denote a set of parameters of the MM+ approach by taking into account the contribution of aromatic resonance as MM+(a). To elucidate the role of aromatic resonance in structural stabilization, the MM+(a) simulations were accompanied by structural optimization without taking into account the aromatic resonance using a set of MM+ parameters denoted by MM+(b). The MM+(b) stabilization energies of n × n flakes and n × 5 nanoribbons are presented in figure 17(f). Since the aromatic resonance almost completely stabilizes the planar topology of 2D biphenylene lattices, the n × n and n × 5 MM+(a) curves have almost constant stabilization energy values which are very close to 0, so they are not displayed in figure 17(f).
Like PM3 stabilization energies, the n × n and n × 5 MM+(b) ones also demonstrate 1/x patterns with much smaller saturation limits of −0.002 and −0.001 kcal/mol/C atom. Once again, these values are much smaller of both vdW energies [112] and the accuracy of the MM+ approach [66,130]. The curvature radii of the largest curved MM+(b) C 726 H 66  The n × 5 MM+(b) nanoribbons reveal a distinctive curved shape ( figure SI1-14). In particular, the longest n × 5 C 630 H 94 nanoribbon of 91.862 Å length has R n×4 = 196.2 Å curvature radius with the largest departure from xy plain equal to 10.337 Å. As one can see, the curvature radius of n × 5 MM+(b) nanoribbons is significantly larger than the PM3 one of R n×4 = 85.139 Å (see above).
Taking into account the aromatic resonance contribution by applying of MM+(a) set of parameters leads to localization of almost perfect planar structures of n × n flakes with maximum departure of carbon atoms from the plane by 0.532 Å, and R.M.S. of D 2h symmetry group of 0.010 (R.M.S. of C 2v symmetry group is 0.000) for the largest curved C 726 H 66 11 × 11 flake with 47.363 × 38.784 Å dimensions. The 7 × 7 and 9 × 9 MM+(a) curved and planar flakes revealed no imaginary frequencies in Hessian calculations.
The MM+(a) set of parameters also leads to localization of almost planar n × 5 nanoribbons with stabilization energies very close to 0, which is much less the accuracy of MM+ approach (not presented in figure 17(f) because of almost perfect coincidence of the curve with X axis). The maximum departure of carbon atoms of the longest (91.316 Å) 20 × 4 nanoribbon from parent planar structure is 1.276 Å with curvature radius equal to 1633.0 Å and D 2h R.M.S. equal to 0.274. It is necessary to note, that the maximum departure of carbon atoms for 15 × 5 nanoribbon is equal to just 0.054 Å with D 2h R.M.S. equal to 0.012.
In table 2 structural and energetic characteristics of n × n, n = 3-9, k × 5, k = 7-15, and 5 × l, l = 7-15 2D biphenylene fragments of finite dimensions calculated using ab initio B3LYP/6-31G * approach are presented. Following the results of PM3 simulations in which the curvature radius was estimated to be 85Å (see above), an extended ab initio B3LYP/6-31G * geometry search was performed starting from curvature radius of 30Å. It is necessary to note that independently of the starting geometry, geometry relaxation without symmetrical restrictions returned the same type of slightly corrugated structures. To reveal the stabilization bending energies and symmetry effects, perfectly planar conformers were calculated as well using the same ab initio B3LYP/6-31G * approach. Table 2. Size (in unit cells and Å), chemical formula, stabilization energies (kcal/mol/C atom), departure from plane (Å), D 2h R.M.S. and nature of planar and curved conformers of n × n, n = 3-9, k × 5, k = 7-15, and 5 × l, l = 7-15 2D biphenylene finite fragments calculated using ab initio B3LYP/6-31G * method a . All 3 types of fragments revealed the bending energies of 0.0-9 × 10 −3 kcal/mol/C atom, which is much lower the accuracy of the DFT B3LYP approach [68,69,[131][132][133][134][135], for which the dissociation energies for the G2 set are determined with an accuracy of 3.5 kcal/mol. Calculated maximum of absolute departure of carbon atoms from the structural plane of 2D biphenylene fragments reaches 0.165 Å for one of 5 × 7 C 210 H 34 conformers of 21.078Å × 25.079Å dimension with D 2h symmetry R.M.S. 0.068 for 9 × 5 C 270 H 46 nanoribbon, which provide unequivocal proof of visible distortion of parent planar D 2h symmetry of all fragments except 3 × 3 C 54 H 18 flake, for which internal total mechanical stress is not strong enough and cannot overcome π-electronic resonance making the lattice perfectly planar. It is necessary to note, that DFT B3LYP accuracy in the determination of interatomic distances is well below 0.001 Å (see references above).
Anizotropic nature of 2D biphenylene lattice leads to visibly different departures from planar structures for k × 5 and 5 × l nanoribbons in perpendicular directions (table 2) equal to 0.070 and 0.165 Å and D 2h R.M.S. deviation of 0.068 and 0.043, respectively. Because of negligibly small stabilization energies of all considered fragments all of which are characterized by different structural shifts in perpendicular directions, one can conclude that 2D biphenylene lattice itself is characterized by very shallow potential energy surface in the vicinity of global minimum, which is close to perfectly planar structure. Since the basic structural units of 2D biphenylene, namely cyclobutadiene and cyclooctotetraene are tensed correlated antiaromatic structures, the entire 2D crystalline lattice of biphenylene could be regarded as a strongly correlated structure [136] which is characterized by both dynamical [137] and non-dynamical (see, for example, [138]) contributions.
The geometry optimization procedures are based on the search of zero gradient on total potential energy surfaces [139]. Rigorously, the outcome of geometry optimization should be accompanied by Hessian calculations at the extremum point to prove whether or not it is a local or global minimum rather than saddle point on complete potential energy surface. In the case of absence of imaginary modes in the calculated vibration spectrum, one can conclude that the localized geometry is a minimum, otherwise, it is proved that it is a saddle point. In general, for extended clusters (hundreds and thousands of atoms) the Hessian calculations are incredibly complex and unreliable making impossible to reveal the true nature of localized extreme points on potential energy surfaces.
For most cases, the optimization trajectories follow the minimum energy path on potential energy surfaces which couple the minimization of the gradient with a smooth decrease of the total energy in the vicinity of an extremal point ( figure SI1-15), which can be either a minimum or a saddle point. Let us denote the first type of trajectories as 'no conclusive answer' (N.C.A.) trajectories since one should run additional Hessian calculations to reveal the true nature of the localized structure. For all N.C.A. trajectories, which couple the local extremum of potential energy surface with gradient minimum, the starting point has higher energy in respect to final localized point. In particular, even the localized extremal point is a saddle point, the N.C.A. trajectory returns a decline of energy in respect to the optimization steps. In general, no conclusive deductions about the nature of the localized point can be made based on just an analysis of the N.C.A.-type optimization trajectory.
Another special outcome of optimization procedure couples the gradient minimization with a raise in total energy ( figure SI1-15). By definition, a trajectory along the raising energy pathway may only lead from a lower energy point on the potential energy surface to a saddle point with higher energy, or global or local maximum. Such sort of trajectories could be denoted as TS trajectories since they may lead to transition states of different orders with higher energies in respect to the starting point. In general, the TS-type optimization trajectory is a clear and unequivocal indication of non-ground state nature of any localized extremum point on potential energy surface as an outcome of optimization procedure.
Step-by-step extension of biphenylene 2D l × m clusters leads to decreasing the influence of cluster boundaries on the flake lattices improving the quality of initial structural guess. Based on TCT analysis (see above) and negligibly small energy differences between planar and non-planar conformers (table 2), one can expect that perfectly planar flakes may be a transition state on potential energy surfaces among multiple almost energy degenerated conformers of undetermined nature. The ab initio B3LYP 6-31G * optimization trajectories are presented in tables 2 and SI1-4. One can see that for extended perfectly planar clusters starting from n × n, n 7, k × 5, k 11 and 5 × l, l 9 the planar conformers were localized within typical TS-type optimization trajectories, which might indicate their transition state character. Even the lattices of some curved conformers ( Coexistence of multiple degenerated or close to degeneracy (energy differences 0-10 −3 kcal/mol/C atom, table 2) correlated transition states, local/global minima or maxima of the same symmetry and very close structural parameters of 2D biphenylene conformers contradicts to fundamental quantum adiabatic theorem [140], making impossible to exist stationary quantum state of the lattice. Following the adiabatic theorem and correlation effects within the lattice (see discussion above), the 2D biphenylene can be considered as a strongly-correlated non-adiabatic 2D crystal. Nevertheless, following the recent consideration [141] of closed quantum systems at final temperatures, in the basis of the instantaneous eigenstates of time-dependent Hamiltonian 2D biphenylene could exist as a correlated mixed quasi-Gibbs quantum state.
The phonon dispersion of perfectly planar 2D biphenylene lattice calculated by the PW PBE PBC method is presented in figure 17(g). An imaginary band with very small absolute amplitude is detected in the vicinity of Γ point. Taking into account the TCT-based structural analysis above one could conclude that the band may be at least partially attributed to a soft mode in the vicinity of global minimum as well as a result of numerical inaccuracy in calculations of phonon dispersion law.
Based on TCT analysis, analysis of possible role of electronic correlations and the effects of non-adiabaticity, it was found that 2D biphenylene could be considered as a strongly correlated non-adiabatic low-dimensional crystal. It was shown, that the aromatic resonance and electronic correlations almost cancel internal structural stress caused by distorted correlated antiaromatic lattice fragments, namely cyclobutadiene and cyclooctotetraene. The structural effects and stabilization binding energy of 2D biphenylene lattice strongly depend upon the method used to simulate the structure. In particular, the PM3 approach [67], which significantly overestimates the force constants, predicts the formation of biphenylene rings with ∼85 Å radius. It is necessary to note that even at PM3 level of theory the stabilization bending energy is smaller than the vdW energy making 2D biphenylene lattice possible to be synthesized atop of some kind of support. Using model potential simulations with and without taking into account the aromatic resonance it was found that aromatic resonance almost completely compensates accumulated structural stress. ab initio B3LYP/6-31G * calculations of extended k × l, k, l = 3-15 also demonstrate negligibly small stabilization bending energies with soft modes in the vicinity of the global minimum. Based on adiabatic theorem consideration, it was shown that free-standing low-dimensional biphenylene lattice can exist at final temperatures as a quasi-Gibbs mixed quantum state. Experimental synthesis of 2D biphenylene by bottom-up approach [28] on solid-state support is in perfect agreement with theoretical results obtained by TCT analysis accompanied by model potential and electronic structure simulations.

TCT analysis of topological and quantum stability of low-dimensional lattices
Following TCT treatment coupled with numerical simulations, a flowchart (figure 18) for TCT structural analysis of low-dimensional lattices might be proposed. At the first step, after a theoretical design of low-dimensional crystalline lattice, one should verify whether or not relaxed free-standing structural units can perfectly fill in the 2D crystalline space. If the fragments perfectly fill in the lattice, PBC calculations of atomic and electronic structure and phonon dispersion can be performed and analyzed. For this particular case, the phonon calculations can be used to prove the stability of proposed low-dimensional materials. However, it is also suggested to carry out ab-initio molecular dynamic simulations to exclude false-predicted nanocrystals [142].
If the proposed lattice does not satisfy TCT requirements, rigorous study of the lattice stability should be performed. In this case, at the second step, conservation of TI conditions should be tested using, for example, electronic structure calculations of finite clusters. If the structural mismatch of the structural units leads to just multiplication of translation vectors in one or both directions, the type and characteristics of regular lattice deformations in a form of a structure wave should be determined (wavelength and amplitude, for example). For regular periodic patterns, the quantum stability of the lattice should be considered and stabilization bending energy should be compared with the energy of van-der-Waals interactions or external pressure to determine the ways to stabilize the proposed lattice by some sort of supports, external pressure, or formation of finite flakes to avoid lattice quantum instability.
Periodicity of the lattice can be eliminated just in one direction with the formation of a nanotube or a nanoroll. In this case, structural characteristics of 1D tubular forms (diameter, chirality, helicity) should be determined.
A complete breakdown of TI conditions can be either achieved in the case of torsion deformation of 2D lattice, by breakdown of the symmetry of the lattice, formation of saddle structures, or in the case of entire decomposition of the closest chemical environment of each single atom. In the case of conservation of the closest environment, small deformed flakes of proposed materials could be somehow synthesized. In this case, the phonon calculations can not be considered as a final and solid proof of structural stability of the lattice because of decomposition of TI conditions. And finally, in the case of very poor initial guess of a lattice, even short-range chemical environment of any constituted atoms may be completely compromised by internal stress with a total breakdown of periodicity. In this case, the phonon calculations can not be considered as final and solid proof of structural stability of the lattice as well and the mechanism of structural instability should be revealed.
The proposed algorithm ( figure 18) can be used to study topological and quantum stability of novel perspective low-dimensional materials.

Conclusions
Short-and long-range orders and deformations of a number of 0D, 1D, and 2D carbon atomic lattices with multiple non-equivalent sublattices were theoretically analyzed and simulated using model potential, semiempirical and DFT approaches. It was found that non-equivalent sublattices may cause uncompensated structural stress and torque with consequent mechanical deformation of the whole lattice. The influence of the key 0D fragments, constituted by 4-, 5-, 6-, 7-, and 8-member rings, on the structure of 2D carbon hexagonal graphene lattice was studied and it was found that all considered cores may cause breakdown of its symmetry and perfectly planar topology.
Based on structural analysis of low-dimensional atomic lattices, a Topology Conservation Theorem was formulated and proved. Assuming that at least one force constant of a low-dimensional crystalline lattice is very small or equal to zero, which is true with high accuracy for one-unit cell thick lattices, it was found that in contrast to conventional 3D crystals, the lack of perfect filling of 2D active crystalline space may lead to either formation of non-planar lattices or even to decomposition of entire regular crystalline structure. It was shown that PBC imply linear translation symmetry restrictions and as a result, causes artificial stabilization of both 1D and 2D low-dimensional lattices prone to violation of perfectly linear or planar topology. In addition, two TCT corollaries were formulated and were shown to be true by comparison with known experimental facts.
All possible types of perfect planar topology violations of low-dimensional lattices with and without breaking of either long-or short-range crystalline orders were considered, namely: (i) spontaneous formation of structure waves or regular helices of different nature with either variable or constant wavelengths with multiplication of translational vectors; (ii) spontaneous formation of 1D nanotubes and rolls; (iii) spontaneous formation of aperiodic saddle structures keeping intact short-range local atomic environment; (iv) complete breakdown of both long-and short-range orders accompanied by a complete breakdown of structural integrity the lattice; (v) possible stabilization of non-planar lattices by any kind of substrates in case of stabilization bending energies smaller the van-der-Waals interactions.
It was shown that the breakdown of TI causes dramatic expansion of complete active space of normal coordinates X for which the TI coordinates correspond to a low-dimensional subspace x ∈ X. Constrained geometry optimization of the lattice in restricted TI coordinates subspace x may lead to localization of ordinary point on complete Hilbert space of normal coordinates as an extremum with consequent erroneous calculations of phonon dispersion law at this point. It was found that for low-dimensional crystalline lattices with multiple non-equivalent sublattices the absence of imaginary modes in theoretical phonon spectra cannot be used as a solid and final proof of their either stability or metastability and should be combined with TCT analysis and robust consideration of finite clusters.
In particular, it was found that notorious 2D phC(m,n) (m;n = 0;1) lattices are prone to form superperiodic structure waves, screws or helices with variable or constant wavelength, 1D nanotubes or nanorolls, and aperiodic atomic lattices as well, with perfectly planar regular 2D lattices as transition states between symmetrically equivalent bent global minima. It was found that at all levels of theory the stabilization bending energies of phC(m,n) lattices are of the same order of magnitude or even larger than the van-der-Waals interactions, which makes it impossible to stabilize planar conformers on top of any kind of supports at finite temperatures.
The standing structure waves for finite-size flakes can be described by homogeneous D'Alembert second-order differential equations for which the atomic positions are well-determined. In contrast to it, infinite free-standing structure waves are described by inhomogeneous D'Alembert equation which has a solution of an infinite superposition of plane waves with non-zero external force, for which the positions of all atoms in the lattice are completely undetermined. It was shown that finite wave conformers can exist and can be synthesized, whereas infinite free-standing wave lattices are prone to lattice quantum instability effects. In particular, the lattice quantum instability should lead to the breakdown of wave periodicity by instantaneous structural transformations.
In contrast to previous publications and in consistency with TCT analysis and finite cluster considerations, the phonon dispersion of planar 2D phC(0,1) lattice is characterized by one imaginary mode. The one-wave-length supercell PW PBC PBE electronic structure calculations of both wave and planar conformers of phC(0,1) unequivocally confirm the global minimum and transition state nature of wave and planar conformers, respectively. The phonon dispersion laws of phC(1,0) and phC (1,1) lattices also revealed imaginary modes in consistency with TCT and finite cluster behavior. Despite the screw behavior of elongated clusters and violation of TCT mandatory requirements, the phonon dispersion law of 2D phC(0,0) R 5,7 Haeckelite does not reveal any imaginary modes, which was interpreted as a result of the breakdown of lattice TI and phonon calculations performed at constrained minimum in TI coordinates subspace of normal coordinates. It was speculated that for free-standing 1D phC(0,0) nanoribbons which follow screw pattern, the lattices should suffer quantum instability.
It was found that at all levels of theory the Hessian calculations of finite phC(m,n) flakes require incredibly high accuracy of minimization of the lattices, with the gradient stopping criterium better than 10 −5 kcal/mol/Å. Nevertheless, some obvious mistakes in calculations of vibration spectra made by different codes for extended low-dimensional clusters were detected.
It was shown that NVT AIMD PBC algorithm artificially confines periodic lattices inside MD simulation box of constant dimensions which prevents true structural relaxation and thermal equilibration in form of structural contraction or stretching with breakdown of topology and symmetry of low-dimensional crystals. For example, in the case of formation of structural waves one should perform AIMD simulations of 0D finite clusters to test stability of their topology and symmetry.
Using ab initio DFT calculations, the stabilization bending energy of 2D penta-graphene was estimated to be 1.9 kcal/mol/atom which is at least as twice as large as the vdW energy. It was shown that mathematically free-standing penta-graphene cannot form any regular 0D, 2D, or 3D geometrical solids like spherical fullerenes, tori, corkscrew tubes, or triple periodic minimal surfaces. Strong bending force constants cause formation of either 0D small saddle-shaped flakes, 1D nanotube of 11.96/14.27 Å inner/outer diameter, or complex aperiodic 1D nanorolls of 65-70Å external diameters. It was shown that PBC symmetry restrictions artificially stabilize planar 2D penta-graphene, preventing breakdown of planar topology of the lattice. Constrained PBC minimization of planar 2D penta-graphene in translation symmetry subspace leads to artificial localization of a regular point on complete active space natural coordinates potential energy surface as a global minimum with a lack of imaginary modes in phonon dispersion. It was concluded that phonon dispersion cannot be used as final and solid proof of either stability or metastability of 2D penta-graphene lattice.
The TCT analysis of ladderane unequivocally indicated topological instability of the lattice due to dramatic departure of its rectangular structural units from either triangular sp 2 , or tetrahedral sp 3 carbon atom coordinations. It was shown that imposing of translation symmetry leads to artificial stabilization of ladderane keeping intact its rectangular periodic lattice. In contrast to PBC calculations of infinite 2D ladderane, the PW PBE cluster optimization revealed its structural instability with complete loss of periodicity, structural integrity and even local symmetry of carbon atoms making up the initial lattice. It was shown that complete active space of ladderane normal coordinates of infinite dimensionality includes translation symmetry subspace, in which constrained minimization of the structure returns artificially stable geometry. The phonon dispersion law obtained at constrained minimum returns no imaginary modes because of significant departure of ladderane translation symmetry subspace from complete active space of natural coordinates true minimum, which corresponds to an infinite irregular ensemble of asymmetric carbon clusters of different nature. It was concluded that the phonon spectrum of ladderane cannot be used as final and solid proof of its structural stability/metastability. Using TCT analysis followed by model potential simulations and electronic structure calculations of 2D biphenylene it was shown that the lattice does not satisfy the mandatory requirements of the theorem. It was found that stabilization of planar 2D biphenylene topology is caused by aromatic resonance and correlation effects which almost completely compensate internal mechanical stress caused by structural inconsistency of the fragments. Based on simulations of atomic and electronic structure it was found that stabilization energy of 2D biphenylene in order of magnitude is lower than the vdW interactions which allows the lattices to be synthesized at any kind of solid support. It was shown that free-standing biphenylene is a correlated non-adiabatic low-dimensional crystal that can exist at final temperatures as a quasi-Gibbs mixed quantum state.
Following the results of theoretical consideration and numerical simulations of 2D carbon lattices it was shown that TCT imposes rigorous mandatory limitations on topology, symmetry and structure of low-dimensional lattices with multiple sublattices. Based on TCT theorem, a flowchart of TCT structural analysis was proposed. The algorithm can be used to develop advanced methods to analyze structure and topology of low-dimensional crystalline lattices.
1        [1] for 22 and 33 supercells calculated using PHONON code [3][4][5]. No imaginary frequencies were detected. Figure SI1-6b. Phonon dispersion for phagraphene [2] modeled by Tersoff potential. No imaginary frequencies were detected. The spectrum was calculated using GULP [6,7] code.     Penta-graphene flakes, rings, tubes and rolls calculated using PM3 approach. a) Comparative side view of penta-graphene × 8 ( = 5 − 8) flakes. 5 × 8 C360 flake is presented in yellow, 6 × 8 C396 flake is presented in red, 7 × 8 C420 flake is presented in green and locked 8 × 8 C430 flake is presented in blue. b) Two side views of 5 × 8 C360 flake. c) Two side views of 6 × 8 C396 flake. d) Two side views of 7 × 8 C420 flake. e) Two side views of 8 × 8 C430 flake. f) and g) Two side views of unlocked and locked of b-oriented ring. h) and i) Two side views of unlocked and locked of b-oriented nanotube fragment of 5 unit cell length. j) A roll formed by rolling of a-oriented nanoribbon. The main axis of the roll is parallel to b-direction. k) and l) Two side views of unlocked and locked of a-oriented ring. m) A roll formed by rolling of b-oriented nanoribbon. The main axis of the roll is parallel to a-direction.

A. 1D h-BN zig-zag narrow nanoribbon case
A hypothetical 1D h-BN zig-zag narrow nanoribbon (h-BN ZNR) of one B3N3 hexagonal fragment width ( Figure SI2-1) may be considered as the simplest case of a low-dimensional lattice with multiple non-equivalent sublattices. Without loss of generality one can introduce a1, a2, and b1, b2 nonequivalent sublattices for N and B basis atoms, respectively, associated with symmetrically nonequivalent aa1, aa2, ab1, ab2 translation vectors oriented along X direction. In particular, both external atomic rows have coordination numbers equal to 2 (or, in the case of hydrogen-substituted external atoms, 3, with symmetrically different environment), with different neighborhood, namely a1 N atoms have 2 b2 B neighbors while b1 B atoms have 2 a2 N neighbors. The coordination numbers of a2 N and b2 B atoms are equal to 3 with 2 b1 and 1 b2 boron atoms and 2 a1 and 1 a2 nitrogen neighbors, respectively. For simplicity all interatomic bond lengths can be considered equal to RNB. Since the basis atoms, the boundary conditions, coordination numbers and types of environments of a1/a2 and b1/b2 N and B sublattices are nonequivalent, one can write: The N a1 and B b1 sublattices with coordination numbers equal to 2 possess different force constants Q1, Q2 of vibrations along X direction because of different symmetry and environment of the edge N and B atoms even a2-b1-a2 and b2-a1-b2 angles and RNB length of all N-B bonds keep initial parameters of perfect 2D h-BN. The force constant Q along the X direction and RNB length of B-N bonds for a2 and b2 sublattices with coordination numbers equal to 3 keep the initial Qi value of perfect 2D h-BN lattice. So, for aa2 and ab2 translation vectors and the force constants the following relationships can be written: sign. In respect of the origin of the coordinate system (red dot in the center Figure SI2-2), the coordinates of atoms 1 and 2 are: In harmonic approximation the X and Y components of a1, a2 associated energies and forces are: This conditions can be satisfied only all components are equal to 0, so, for pristine hexagonal lattices of graphene and h-BN, the X and Y components of forces mutually compensate each other and the total stress of the lattice is equal to 0. For the torques acting on atoms 1 and 2 in respect to the origin of the coordinate system one can write:

Estimation of internal pressures and temperatures of phC(0,1) structural waves
For phC(0,1) nA✕1B lattice two significantly different structural waves were localized. The first one, of 5 unit cell wavelength and 1 unit cell width is characterized by 33.8Å wavelength, 6.9Å width and 9.2Å 2 ⁄ = 4.6Å amplitude (see Figure S4). The volume of a box which contains one structure wave unit is: 7 = 33.8 × 9.2 × 6.9 × 10 −30 3 = 2.1 × 10 − 27 3 Each single structure wave box contains 137 carbon atoms, so each single atom occupies = 2.1 × 10 −27 3 137 ⁄ = 1.5 × 10 −29 3 . Averaged carbon-carbon bond length is 1.5Å, so a carbon atom occupies 1. Under normal conditions, each single atom of an ideal gas occupies 3.7 × 10 −26 3 . Since the atomic positions of an infinite structural wave are completely uncertain due to quantum effects, one can speculate that atoms are evenly distrubutted in the structural box determined by the effective wave width, wavelength and amplitude. The effective internal pressure of acting on each atom of the unconstrained and undisturbed infinite structural wave can be calculated as a quotient of the atomic volumes under normal conditions and in the structure wave: 3.7 × 10 −26 3 1. The second structural wave of 18 unit cell wavelength is localized for phC(0,1) nA4B nanoribbons and nAnB (see above, Figure S5). It has 115.8 Å wavelength, the smallest width of with averaged atomic volume of 6.5 × 10 −29 3 which corresponds to 569.2 atm pressure and averaged carbon-carbon distance of 4.0Å.
Once again, using Ideal gas law for 18A4B structural wave the effective temperature is: 8

Estimation of bending energy of penta-graphene clusters with taking into account of boundary effects
For penta-graphene 3 × 3 C72H28 and 4 × 4 C120H36 clusters the boundary effects can be taken into account supposing that two distinctively different types of atoms at the edge of the clusters with reduced and preserved coordination numbers contribute differently to stabilization energy. Using linear approximation one can simply write two independent linear equations where Es and Eb are stabilization energies for small C72H28 and extended C120H36 clusters, respectively, es and ee are strain energies for perfectly coordinated and edge atoms, respectively, 1 = 26 and 2 = 36 are the numbers of edged carbon atoms for C72H28 and C120H36 clusters, respectively, and 1 = 46 and 21 = 84, for C72H28 and C120H36 are the number of carbon atoms with preserved coordination numbers. Solution of the system returns B3LYP/6-31G * strain energy for non-edged atoms = 1.879 ⁄ / , which is way beyond vdW binding energy of 9.6×10 −1 − 9.6×10 −2 kcal/mol/atom, [3].