Mechanical Properties of Cubene Crystals

The fullerene family, whose most popular members are the spherical C60 and C70 molecules, has recently added a new member, the cube-shaped carbon molecule C8 called a cubene. A molecular crystal based on fullerenes is called fullerite. In this work, based on relaxational molecular dynamics, two fullerites based on cubenes are described for the first time, one of which belongs to the cubic system, and the other to the triclinic system. Potential energy per atom, elastic constants, and mechanical stress components are calculated as functions of lattice strain. It has been established that the cubic cubene crystal is metastable, while the triclinic crystal is presumably the crystalline phase in the ground state (the potential energies per atom for these two structures are −0.0452 and −0.0480 eV, respectively).The cubic phase has a lower density than the monoclinic one (volumes per cubene are 101 and 97.7 Å3). The elastic constants for the monoclinic phase are approximately 4% higher than those for the cubic phase. The presented results are the first step in studying the physical and mechanical properties of C8 fullerite, which may have potential for hydrogen storage and other applications. In the future, the influence of temperature on the properties of cubenes will be analyzed.


Introduction
The existence of large carbon molecules was predicted and described in some early papers, but a paper [1] published in 1985 caused an avalanche of publications on fullerenes. The smallest spherical fullerene is C 20 ; the most common molecules are C 60 and C 70 , but fullerenes with a larger number of carbon atoms are easily obtained [2,3]. Probably the smallest possible carbon cluster is a cubic C 8 molecule called cubene [4].
Cubenes C 8 and cubanes C 8 H 8 [5,6] are less studied than C 60 and C 70 fullerenes. The cubene molecule is highly strained, since a bond angle of 90 • is energetically unfavorable; however, as follows from the ab initio study [7], the molecule should be stable because its energy corresponds to a deep minimum. Structural, electronic, and optical properties of cubane C 8 H 8 and cubanoids (cubane-like molecules with hydrogen substituted by other atoms) were analyzed using a DFT approach [8].
Derivatives of cubane such as octanitrocubane and heptanitrocubane are excellent candidates for storage of large amounts of energy [9]. Thermochemical calculations show that 4145 kJ/mol is released during decomposition of the first, and 4346 kJ/mol of the second. Cuban can be used to obtain various bioisosteres, i.e., compounds obtained by replacing a group of atoms with a cubane, while maintaining the biological activity of the original compound. Thus, new substances with improved properties compared to the parent compound for pharmaceutical and agrochemical applications can be created [10]. Cubanes are also promising for hydrogen storage technologies [11][12][13].
Covalently bonded arrays of cubanes were shown to be possible [14]. Chemical modification of cubanes and cubane oligomers was analyzed in the work [15]. The electronic structure of C 8 and B 4 N 4 cubic systems was analyzed [16]. It was shown that the functionalization of a chain of covalently bonded cubanes leads to a significant change in its conductivity [17,18]. This means that the electronic properties of cubane-based oligomers can be tuned by functionalization for applications in nanoelectronics.
Molecule crystals composed of fullerenes are called fullerites; they possess many promising physical, mechanical, and chemical properties for applications. Polymerization of fullerites under high pressure and temperature produces extremely hard materials with bulk modulus from 454 GPa to 644 GPa [19]. Fullerite demonstrates an anomalously slow thermal relaxation due to the excitation of internal vibrational degrees of freedom [20].
The crystal structure of fullerites is of great importance, since the structure determines the properties of crystals. Fullerite C 60 can be packed into a simple cubic or face-centered cubic (fcc) lattice [44,45]. In both lattices, a primitive translational cell contains one fullerene molecule. In the theoretical work [46], it was shown that, for C 60 , the fcc structure is more stable than hexagonal close-packing (hcp), but for C 70 hcp has lower energy than fcc lattice. The fcc structure of C 60 at 0 K has an orthorhombic distortion due to ordering of fullerenes. At higher temperatures, there is a transition to an orientationally disordered fcc lattice [46]. As can be seen, many studies have been carried out to predict and analyze the structure of fullerites C 60 and C 70 , but we have not been able to find any information about the structure for fullerite C 8 except for the study [47] where cubic cubane C 8 H 8 was described. In the present study, we set ourselves the task of finding the equilibrium phases of a cubene crystal and their mechanical properties.

Materials and Methods
Molecular dynamics is used to simulate the structure and mechanical properties of cubene crystals. Carbon-carbon interactions are described by the harmonic potential, while non-valence interactions are described by the Lennard-Jones potential. Each cubene interacts with the nearest, second and third neighboring cubenes. Below is a detailed description of the model.
A model of cubene molecule is shown in Figure 1. Each vertex of a cube with side ρ 0 contains a carbon atom. The covalent bonding of atoms in cubenes makes them relatively rigid. Cubenes interact with each other through much weaker van der Waals bonds.
In the present work, the equilibrium phases of a cubene crystal are described, and, in this case, cubene deformation plays an insignificant role. Therefore, covalent bonds are described by a simple harmonic potential. It is assumed that each carbon atom interacts through covalent bonds with the three nearest atoms and with the most distant atom, for example, atom 1 interacts with atoms 2, 4, 5, and 7 (see Figure 1). The interaction between the farthest neighbors is introduced to ensure the stability of the cubic shape of the molecule. Numbering of atoms in a cubene is shown. It is assumed that each carbon atom interacts through covalent bonds with the three nearest atoms and with the farthest atom; for example, atom 1 interacts with the atoms 2, 4, 5 and 7. Covalent bonds are described by the harmonic potential of stiffness K and equilibrium length ρ 0 for the nearest neighbors and √ 3ρ 0 for the farthest neighbors.
The equilibrium phases of cubene crystal are found using perturbative relaxational molecular dynamics. In this approach, the positions of the atoms are first perturbed to drive the system away from a local potential minimum, and then the structure is relaxed to reach a deeper energy minimum. Initially, cubenes with random orientation are arranged in the fcc lattice with periodic boundary conditions. The simulations are carried out for a computational cell that includes 2 × 2 × 2 primitive translational cells of the fcc lattice. The total number of cubenes in the computational cell is M = 8, and the total number of carbon atoms is 64. The structure relaxation is performed with respect to the coordinates of the carbon atoms and with respect to the vectors specifying the periodic cell of the structure in order to achieve the state with zero forces acting on the atoms and zero external stresses applied to the computational cell.
Our model is built in the spirit of the work [20], where the deformation of valence bonds is described by the harmonic potential, and intermolecular interactions by the Lennard-Jones potential. Following this work, the potential energy of the cubene crystal described by the periodic computational cell with M cubenes can be written as follows: The first term on the right-hand side of Equation (1) describes the valence interactions between nearest atoms in the cubenes; there are 12 nearest bonds in each cubene. The harmonic potential is used: where K = 31.77 eV/Å 2 is the bond stiffness [48], |r n | is the distance between two nearest neighbors in a cubene molecule and ρ 0 = 1.418 Å is the equilibrium bond length. From Figure 1, it can be seen that the twelve bonds between nearest neighbors connect the following pairs of atoms: 1-2, 2-3, 3-4, 4-1, 1-5, 2-6, 3-7, 4-8, 5-6, 6-7, 7-8 and 8-5. The second term on the right-hand side of Equation (1) describes the valence interactions between farthest atoms in the cubenes; there are four longest bonds in each cubene; they are described by the harmonic potential: where K is the bond stiffness the same as for the bonds between nearest atoms, |r k | is the distance between two farthest neighbors in a cubene molecule and √ 3ρ 0 = 2.456 Å is the equilibrium bond length. From Figure 1, it can be seen that the four bonds between farthest neighbors connect the following pairs of atoms: 1-7, 2-8, 3-5 and 4-6. The third term on the right-hand side of Equation (1) describes the van der Waals interactions between carbon atoms of different cubenes. The computational cell includes 8M atoms numbered by the index i = 1, 2, ..., 8M. The i-th atom has the radius vector r i = (r i,x , r i,y , r i,z ), and R i,j = r j − r i . The van der Waals interactions are described by the Lennard-Jones potential where = 0.002757 eV and σ = 3.393 Å [49]. The van der Waals interactions are taken into account with the cut-off radius so that each cubene interacts with 12 nearest, 6 next-nearest, and 24 next-next-nearest cubenes. Van der Waals interactions between atoms belonging to the same cubene are not taken into account. A comment should be made on the choice of interatomic potentials used in this work. To date, a number of interatomic potentials have been developed for carbon materials, the most popular of which are REBO-I [50], REBO-II [51], C-EDIP [52], Tersoff [53], LCBOP-I [54], ReaxFFC-2013 [55] and AIREBO [56] potentials and recently developed machinelearning potential GAP-20 [57]. Each potential has been designed for specific purposes. In this study, we do not consider high pressures or temperatures, so the potential can be greatly simplified to speed up the simulation. The van der Waals bonds describing the interaction between cubenes are two orders of magnitude weaker than covalent bonds within a cubene. Therefore, for our study for covalent bonds, the simplest harmonic potential is sufficient.
A home-made code written in the C++ algorithmic language is developed for simulations. Components of the stress tensor σ ij are calculated as described, e.g., in [58]. The gradient method is used to minimize the potential energy of the crystal Equation (1) with respect to the coordinates of carbon atoms and with respect to the translation vectors of the computational cell. The relaxation is complete when the maximal force acting on atoms does not exceed 10 −9 eV/Å and the value of the maximal stress component does not exceed 10 −6 eV/Å 3 .
For the relaxed structure, the potential energy per atom and the elastic constants C ij are calculated as described, for example, in [58]. By doing so, the elastic constants in Hooke's law  are determined. Note that, for any crystal, C ij = C ji and, consequently, the maximum number of elastic constants is 21, although there are fewer of them for crystals with high symmetry. We also calculate the potential energy, stress components, and elastic constants of the cubene crystal as the functions of homogeneous strain. After applying the homogeneous deformation, the structure is relaxed with respect to atomic displacements in order to achieve practically zero forces acting on the atoms. Stress components and elastic constants are calculated for relaxed structure. Briefly, the research flow can be described as follows. First, equilibrium phases are found, then the stability of the phases with respect to uniform deformation is analyzed, and then properties of the phases, such as energy, stress-strain curves, and elastic constants, are calculated.

Results
Let us present the numerical results obtained for cubene crystals. First, two equilibrium phases will be described. Then, the stability of the phases with respect to uniform elastic deformation will be demonstrated. Finally, elastic constants and stress-strain curves will be presented for both equilibrium phases.

Two Equilibrium Phases
The structure relaxation was carried out for 50 initial configurations with randomly oriented cubenes having centers of gravity at the sites of the fcc lattice with the distance between nearest neighbors equal to 5.3 Å. In many cases, metastable structures with a relatively high potential energy have been obtained. In some cases, two equilibrium phases with a relatively low potential energy were found. One of them has cubic symmetry, and the other belongs to the triclinic system. These two phases are described below.
The structure of the crystal will be defined by the translation vectors of the primitive cell, a, b, c, and coordinates of the carbon atoms in the cell, ξ i . Then, the radius-vector of any atom in the crystal can be expressed as follows: where α, β, and γ are arbitrary integers.

Metastable Cubic Phase
A primitive translation cell of cubic phase includes two cubenes, see Figure 2. Translation vectors of a primitive cell for the cubic phase are where δ = 7.  Vectors of the primitive translation cell, a, b, and c, are given by Equation (7). A primitive translational cell includes two cubenes differently oriented in space.

Stable Triclinic Phase
A primitive translation cell of triclinic phase includes one cubene, see The volume per cubene is 97.70 Å 3 , which is less than that of the cubic phase and, therefore, the triclinic phase is denser.   (10). A primitive translational cell includes one cubene; therefore, all cubenes have the same orientation in space.

Energy of Cubic and Triclinic Cubene Crystals
Potential energy per atom for equilibrium cubic cubene is −0.045185 eV and, for equilibrium triclinic cubene, it is lower and equal to −0.047987 eV. The triclinic phase has the lowest energy among all structures obtained by relaxation of randomly oriented cubenes. We assume that the triclinic cubene crystal has the lowest possible energy, but we cannot provide a rigorous proof of this fact.
It is important to demonstrate that the crystals are stable with respect to application of elastic strain. The potential energy per atom was calculated for various values of the strain component ε ij , while all other strain components were equal to zero. The results are plotted in Figure 4 for (a) cubic and (b) triclinic crystals. Note that the application of any strain component, positive or negative, leads to an increase in the potential energy of the crystal. This means that cubic and triclinic crystals are structurally stable up to ±5% strain. . For a cubic crystal, the curves for the normal strain components ε xx , ε yy , and ε zz in Figure 4a coincide, the same is true for the shear strains ε xy , ε xz , and ε yz . This fact reflects the lattice symmetry.
For the low-symmetry triclinic phase, the curves in Figure 4b do not match, but the difference is small for normal and shear strains. It can be concluded that the triclinic phase is weakly anisotropic.

Elastic Constants for the Cubene Crystals
The stress components are related to the strain components according to Hooke's law Equation (5), which for the cubic crystal assumes the form (elastic constants are given in GPa)  It is well known that isotropic elastic media are characterized by two independent elastic constants, while crystals with cubic symmetry generally have three independent elastic constants [59][60][61][62]. There are only two independent elastic constants in Hooke's law Equation (13), since pairwise interatomic potentials are used, and the Cauchy relation is satisfied, which imposes an additional restriction on elastic constants [63] .
For the triclinic cubene crystal, the Hooke's law reads (elastic constants are given in GPa)  In Equation (14), due to the Cauchy relations, one has C 44 = C 12 , C 55 = C 13 , and C 66 = C 23 . Analysis of the elastic constants confirms that the triclinic cubene crystal is weakly anisotropic. The constants C 11 , C 22 , and C 33 differ by no more than 1.6%. The constants C 44 , C 55 , and C 66 differ by no more than 4.2%.
A comparison of the elastic constants of the cubic and triclinic cubenes indicates that the latter is stiffer. Indeed, in Equation (13), C 11 = C 22 = C 33 = 464.5 GPa, and, in Equation (14), these elastic constants are higher. The same applies to the shear moduli, which, for a cubic crystal are C 44 = C 55 = C 66 = 177.1 GPa, which is less than for a triclinic crystal. The higher stiffness of the triclinic crystal is due to its greater density.

Stress-Strain Curves
The stress-strain curves for cubic and triclinic cubene crystals are presented in Figures 5 and 6, respectively. For the abscissa, the only nonzero strain component is used, and the other strain components are zero. The normal stress components are shown by black (σ xx ), red (σ yy ), and blue (σ zz ) lines. The shear stresses are plotted by yellow (σ xy ), cyan (σ xz ), and gray (σ yz ) lines. In some cases, dashed lines are used to show overlapping curves.
As can be seen from Figures 5 and 6, when normal strain is applied, the shear stresses are equal to zero in cubic crystal, but, in the triclinic crystal, they are not.
For the normal stress components, as the functions of the normal strain components [panels (a) to (c) in Figures 5 and 6] within the studied strain range ±5% demonstrate, there is a noticeable deviation from the linear dependence. On the other hand, the stress-strain curves σ ij (ε ij ) for i = j show a linear relationship in this strain range (see panels (d) to (f) in Figures 5 and 6).
The elastic constants C 11 , C 22 , and C 33 as the functions of strain are plotted in Figures 7 and 8 for cubic and triclinic cubene crystals, respectively. In the studied strain range ±5%, these elastic constants depend almost linearly on the normal strain components and weakly depend on the shear strain components.
The elastic constants C 44 , C 55 , and C 66 as the functions of strain are plotted in Figures 9 and 10 for cubic and triclinic cubene crystals, respectively. In the studied strain range ±5%, these elastic constants depend almost linearly on the normal strain components and weakly depend on the shear strain components. Figure 5. Stress-strain curves for cubic cubene crystal. Only one component of the strain tensor is varied, and the rest are equal to zero. The normal stress components are shown by black (σ xx ), red (σ yy ), and blue (σ zz ) lines. The shear stresses are plotted by yellow (σ xy ), cyan (σ xz ), and gray (σ yz ) lines. In some cases, dashed lines are used to show overlapping curves. Panels (a-f) correspond to changing ε xx to ε yz , respectively.  . Elastic constants C 11 , C 22 , and C 33 as the functions of strain for cubic cubene crystal (black, red, and blue lines, respectively). Only one component of the strain tensor is varied, the rest are equal to zero. Panels (a-f) correspond to changing ε xx to ε yz , respectively.  Figure 7, but for the triclinic crystal. Panels (a-f) correspond to changing ε xx to ε yz , respectively. Figure 9. Elastic constants C 44 , C 55 , and C 66 as the functions of strain for cubic cubene crystal (yellow, cyan, and gray lines, respectively). Only one component of the strain tensor is varied, the rest are equal to zero. Panels (a-f) correspond to changing ε xx to ε yz , respectively. Figure 10. The same as in Figure 9, but for the triclinic crystal. Panels (a-f) correspond to changing ε xx to ε yz , respectively.

Discussion
Two equilibrium phases of the C 8 molecular crystal were found for the first time based on perturbative relaxation molecular dynamics simulations. The phase of cubic symmetry is metastable with a sublimation energy of −0.045 eV per atom. A primitive translational cell of a cubic cubene contains two cubenes with the volume of 101 Å 3 per one cubene, see Figure 2. A triclinic cubene has a sublimation energy of −0.048 eV per atom, its primitive translational cell includes one cubene, and each cubene occupies a volume of 97.7 Å 3 , see Figure 3. The triclinic phase is most likely the ground state of the cubene crystal.
Stress-strain curves for cubic and triclinic cubenes were calculated in the strain range ±5%, see Figures 5 and 6, respectively. No phase transformations were observed in this strain range; therefore, both phases in this strain range are structurally stable. Components of normal stress in the strain range ±5% show noticeable deviation from the linear Hooke's law.
The elastic constants C 11 , C 22 , and C 33 as the functions of strain for cubic and triclinic cubene crystals are presented in Figures 7 and 8, respectively. Shear moduli C 44 , C 55 , and C 66 depending on strain for cubic and triclinic cubene crystals are shown in Figures 9 and 10, respectively. Elastic constants C 11 , C 22 , and C 33 depend almost linearly on strain, while constants C 44 , C 55 , and C 66 are weakly dependent on strain.
Elastic properties of cubenes can be compared to that of fullerenes. Young's modulus of the fcc single crystal fullerite C 60 at room temperature was measured in [64] to be 20 ± 5 GPa. For C 60 fullerene nanowires, the Young's modulus is 15.1 ± 1.39 [65], while, in bcc tetragonal C 60 nanowire, it is 30 GPa [66]. According to our calculations, stiffness of cubenes is one order of magnitude higher, see Equations (13) and (14). This difference is related to the small size of C 8 cubenes, which greatly increases the number of van der Waals bonds per atom. In addition, note that our calculations are made for zero temperature and, as is well known, the rigidity of molecular crystals decreases very rapidly with increasing temperature [64].
The limitations of our study are as follows. Since the harmonic potential is used to describe covalent bonds, high pressure and the possible transformation of the structure due to polymerization cannot be modeled within our model. The influence of temperature was not taken into account and, consequently, the range of thermal stability of the obtained structures was not established.

Conclusions and Future Work
Two equilibrium phases of C 8 fullerite crystal were discovered by the relaxational molecular dynamics method, one of which has cubic and the other triclinic symmetry, see Figures 2 and 3, respectively. The energy per atom for the triclinic phase is 0.0028 eV lower than that for the cubic phase, and triclinic fullerite is presumably the ground state structure. Both phases are stable with respect to uniform strain ±5%, see Figure 4. The stability of structures with respect to thermal vibrations was not analyzed in this work.
Stress-strain curves for cubic and triclinic fullerite demonstrate noticeable nonlinearity in the range of uniform normal strain ±5%, see panels (a)-(c) in Figures 5 and 6. On the other hand, stress-strain curves for shear strain are almost linear, see panels (d)-(f) in Figures 5 and 6.
The elastic constants C 11 , C 22 , and C 33 as the functions of strain are presented in Figures 7 and 8 for cubic and triclinic cubenes, respectively. Similar results for shear moduli C 44 , C 55 , and C 66 can be seen in Figures 9 and 10 for cubic and triclinic cubenes, respectively. Elastic constants C 11 , C 22 , and C 33 vary almost linearly with strain, while constants C 44 , C 55 , and C 66 are weakly dependent on strain.
The results presented here are the first step in studying the physical properties of cubenes. Future studies will analyze phonon dispersion curves, nonlinear dynamics of cubenes including rotobreathers [31][32][33][34][35], and crowdions [67]. To solve these problems, a more complex interatomic potential should be used [50][51][52][53][54][55][56][57]. The effect of temperature on the mechanical properties and thermal expansion will be analyzed. Another interesting line of future research is the development of continuum models capable of describing the vibrational spectra of cubenes. We hope that our study will stimulate experimental work on the mechanical and physical properties of cubenes. In particular, it would be of interest to analyze the polymerization of cubenes during high-temperature treatment under high pressure, similar to what was done in [68] for C 60   Data Availability Statement: The data of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.