Molecular dynamics simulations of nanosheets of polymeric carbon nitride and comparison with experimental observations

ABSTRACT A computational study of the properties of polymeric carbon nitride using molecular dynamics is presented. The analysis of ideal infinite-extent sheets permits to evaluate the effect of temperature on the network of hydrogen bonds responsible for the linkage of the material. The weakening of this binding mechanism at sufficiently high temperatures, together with the inter-layer interactions characteristic of this type of 2D materials, is shown to determine the conformation properties of polymeric carbon nitride at the nanoscale. The results obtained from the simulation of finite samples in the canonical ensemble at varying temperatures are consistent with those from the characterization of our experimentally synthesized samples. Hydrogen bonding between adjacent polymer ribbons leads this process and is the cause of the typical crumpled structure of this material.


Introduction
Polymeric carbon nitride has attracted large interest as a precursor for the synthesis of super-hard carbon nitride phases. [1][2][3][4][5][6][7][8][9][10] This material finds many applications in the wide field of sensors due to the combination of its chemical sensitivity with its optical and semiconductor properties. [11][12][13] Also, monolayer polymeric carbon nitride forms a new family of wide-gap semiconducting materials, in which the gap can be tuned by small variations in the stress applied to the layers or by the presence of adatoms. [14] The various existing synthesis routes provide predominantly amorphous samples structured at the nanoscale level as graphite-like, [15] nanosheets, [16][17][18][19] nanotubes, [20] or nanopillars. [21][22] They can also be structured in mesoporous conformations that lead to new applications, like the generation of nanostructured metal nitrides. [23] Conformation as single sheets that can be exfoliated from graphitic carbon nitride samples [24] or directly grown [25] permits to increase its efficiency as a catalyst. The detailed arrangement of this material is of most relevance for the specific applications and strongly depends on the ambient conditions of the synthesis [16][17][18][19] and on other parameters, such as the precursor type. [25] Understanding of the detailed supramolecular mechanisms that rule the conformation of this material is thus essential for the optimization of the synthesis procedures and the development of new related materials and applications.
In this work, we perform a study of the conformation properties of polymeric carbon nitride sheets under different ambient conditions using atomistic molecular dynamics (MD) simulations. The analyses address both the general intrinsic properties of ideal infinite-extent 2D polymeric carbon nitride sheets and the evolution in the structural organization of finite size samples; specifically, the routes towards the formation of crumpled nanosheets.
MD [26] is a useful tool in the study of the conformational transitions of very large molecular and supramolecular systems like, for instance, protein folding. [27] It provides not only the final minimum energy structures, but the MD trajectories also afford the detailed pathways and a valuable insight into the underlying mechanisms determining the evolution. Molecular dynamics simulations have been extensively used in the study of the intrinsic properties of 2D materials, particularly graphene [28][29][30][31] and related materials like graphane [32] or multilayer graphene. [33] These systems are especially interesting from a theoretical point of view since, according to the Mermin-Wagner theorem, [34] long-wavelength fluctuations destroy any crystalline order in pure 2D systems. Nevertheless, 2D membranes can exist in a 3D space, where the displacements are not constrained in two dimensions.
Graphene-like materials can be modelled as polymerized or tethered membranes, which have been subject of intense theoretical investigations [35][36][37] and which are characterized by the fact that bonds linking the membrane sites are unbreakable. Even within the physical constraint of self-avoidance, these 2D membranes typically undergo a transition to a crumpled state (or a tubular state in the presence of anisotropy [38] ), whereas a flat state, with strong height fluctuations, is also possible due to the coupling between bending and stretching modes. On the other hand, polymeric carbon nitride sheets are characterized by their binding through a network of relatively weak hydrogen bonds, far from the ideal unbreakable condition, that links the melon linear arrays of the basic heptazine building blocks. Using MD simulations, we find that the morphological properties of the material under study are predominantly determined by the behaviour of the hydrogen bonds under the varying ambient conditions and the inter-planar interactions between the heptazine units. In spite of the limited size of our simulation cells when compared with the actual spatial extent of the investigated features, the results of the computations are consistent with the experimental results.

Experimental
X-ray diffraction measurements X-ray diffraction patterns were obtained by means of a powder diffractometer ENRAF-NONIUS FR590, detector INEL CPS 120 (BRUKER AXS, The Netherlands), with monochromatorquartz discriminator, Debey-Scherrer geometry, using Cu Ka radiation, and glass capillaries for sample mounting. The samples were ground in an agate mortar and sifted. The measurements always lasted for 1 hour, and crystalline silicon was used as a standard.

TEM and SEM characterization
Transmission electron microscopy (TEM) and scanning electron microscopy (SEM) analyses were performed with a JEOL JEM-FS2200 HRP and a JEOL JSM-820 SEM microscope (JEOL, Japan), respectively.

Computational
The MD simulations were performed using the large-scale atomic molecular massively parallel simulator (LAMMPS) [39] using the DREIDING force field. [40] This generic model includes parameters for H and the non-metallic main-group elements, plus a few metals. It specifically takes into account the effect of hydrogen bonds, which are particularly relevant for the material under analysis.
The data used for the MD simulations were built from the 2D ground state geometry of the PM6 Hamiltonian, [41] calculated using MOPAC. [42] The DREIDING force field polymer layers used as initial conditions for the MD computations where relaxed using LAMMPS, both as finite size and periodic infinite-extent sheets. The minimized geometries yielded a D 19.240 A , b D 13.365 A and g D 90.000 . The values of b and g were in excellent agreement with those derived from semiempirical quantum chemistry calculations and previously reported in, [19] with a somewhat larger deviation in the value of a. MD optimized geometries displayed very good planarity, very similar to that of the semi-empirical computations. Two nearly square sheets were generated with 7776 and 31104 atoms and spanning 13 £ 9 and 24 £ 18 unit cells of the 2D polymer, respectively. Their approximate sizes were 159.8 A £ 159.5 A and 319.7 A £ 319.0 A , respectively.
The MD simulations of infinite 2D sheets were performed in the isothermal-isobaric (NPT) ensemble with periodic boundary conditions. The calculations for finite size systems, were started with well-formed one, three or six nearly square polymer sheets, and were performed in the NVT ensemble. The equations of motion used in LAMMPS were those of Shinoda, Shiga, and Mikami [43] and the time integration schemes closely followed those derived in. [44] Thermostating and barostating were implemented using Nos e-Hoover algorithm.

MD simulations
The properties of single sheets of polymeric carbon nitride are first analysed. These are the basic building blocks of three-dimensional arrangements [15,19] and also the desired final products for many applications. These arrangements can potentially be obtained by exfoliation of more complex materials or by direct growth. The study of the response of infinite sheets to changes in temperature and pressure permit to quantify, in particular, the behaviour of the hydrogen bond network as the temperature grows. The results obtained are consistent with subsequent studies on the processes responsible for the morphological organization of the material under treatment at high temperatures.

Infinite sheets
The intrinsic properties of ideal, infinite extent, 2D arrangements were addressed using periodic boundary conditions in the NPT ensemble. Final trajectory equilibrium states for simulations at zero pressure and T D 100 K and T D 700 K are displayed in Figure 3. The results show how a finite temperature induced height fluctuations in the sheet, similarly to the case of graphene [32][33] . The average absolute deviations from the surface level h are shown in Figure 2(a) for the range between 100 K and 1000 K. The results for the two sheets of different size were very similar, with the same approximate linear dependence of h in that range. As opposed to graphene, these height fluctuations have an important effect on the structural integrity of the layer through their impact on the hydrogen bonds network, as it can be appreciated by comparing subplots (a) and (b) with (c) and (d), respectively, in Figure 1. This can be clearly appraised in the dependence of the energy contributed by hydrogen bonds normalized with the total number of heptazine units in the simulation domain shown in Figure 2(b). These results show some dependence with the size of the sample under study. The contribution of the hydrogen bond network to lower the energy of the system decreases as temperature is increased and part of these bonds are broken. Further, the rate of change of this energy with T saturates at higher temperatures due to the finite number of total hydrogen bonds.
The structural integrity of other non-hydrogen-bonded 2D materials is preserved under very large in-plane uniaxial or biaxial stress, resulting in a net change of the potential energy in the membrane. [32] In the polymeric carbon nitride presented herein, simulations conducted even at low temperature show that increasing the applied finite uniaxial or biaxial in-plane stress eventually results in the collapse of the 2D polymer at a finite time that depends on the applied pressure. The atom displacements in the onset of the instabilities resemble those associated to the low frequency lattice modes of the 2D polymer, [19] which may play a key role in the nascent growth of the instabilities.

Finite size samples
The simulations performed on finite size samples show that, in general, the dynamic trajectories leading to their final conformations are driven by two main effects: the previously described weakening of the hydrogen bonds network as temperature increases, plus the plane-to-plane interactions that take place after polymer sections are loosened or released from their hydrogen bond restraints. These two combined effects produce the folding and piling of the material layers in different manners, determining the final arrangement of the material. The detailed evolution depends not only on the ambient conditions, but also on the amount and the characteristics of the initially interacting material.
The conformation dynamics of a finite single sheet of polymeric carbon nitride were analysed from computations performed in the canonical ensemble. The final states of the calculations of 500 psec evolutions of a layer with 13 £ 9 unit cells in a volume of 200 £ 200 £ 200 A 3 at different temperatures are shown in Figure 3. In all cases, equilibration of the structures took place within the first 100 psec.
For the lowest temperatures in the sequence, 500 K and 600 K, the sheet essentially kept its integrity. At 500 K only a very small defect in the 2D structure was found, localized in one of the corners of the sheet, and corresponded to the disordered two-layer piling of a section spanning very few heptazine units. At 600 K, a defect (a bilayer stripe) could be identified. The dynamics of its formation, starting with the folding of one of the corners and the sliding of the feature along the 2D structure, revealed a very interesting balance between the inter-heptazine interaction and its hindrance by the hydrogen bonds at this temperature.
At 700 K, the repetition of the folding and propagation of defect stripes produced the packing of several layers. This disordered and non-homogeneous piling, in turn, induced some curvature in the formed multilayer structures. As the temperature was increased, the folding of the structure and the packing of multiple layers was facilitated by the debilitation of the network of hydrogen bonds. The route to the layer piling then included the break-up of stripes of material and their layer-tolayer interaction. As the hydrogen bond network weakened with the temperature growth, the prominence of inter-planar interaction increased and pilings with larger numbers of layers were produced.
When similar computations were performed on a planar system spanning 26 £ 18 unit cells, this larger sheet seemed to  be more robust against the folding process at 700 K, with morphological modifications that were limited to a somewhat smaller fraction of the total sheet surface than in the previous case. This can be explained by the fact that the initiation of the dynamical evolution eventually leading to the stacking and crumpling of the sheet would be facilitated by the irregularities in the system and the processes would become less likely in the more homogeneous system. At higher temperatures (900 K), the initial breakup of the hydrogen links between 1D chains and their inter-planar interaction led to the formation of a thread of 1D chains stacked in the direction perpendicular to the heptazine planes. Figure 4 shows the results obtained from the equilibration in the NVT ensemble of three 26 £ 18 sheets, which were initially widely separated and with a random disposition. The computations were performed at different temperatures and, in all cases, the duration of the simulation is 500 psec. Equilibrium was reached within the first 200 psec. At the lower temperatures, 500 K and 700 K, the evolution was dominated by the interaction between the well-formed layers, which resulted in their piling in an essentially three-layer system. When the temperature was increased to 900 K, the thermally-driven sheet undulations broke parts of the hydrogen-bond networks, loosening portions of the initial sheets. The subsequent plane-toplane interactions of the fragmented material resulted in the formation of an irregular flake with a curved profile. The number of layers at some locations of this disordered cell was larger than three.
The results derived from the evolution in the canonical ensemble were very similar when the initial conditions consisted of six identical, well separated and placed at random sheets. The time required to reach an equilibrium state in this case was significantly larger. As in the previous case of three sheets in the starting conditions, the evolution at the lower temperatures produced the piling of the initially wellformed layers. The system dynamics, in this regime, were dominated by the layer-to-layer interactions. At 900 K, the thermally induced shaking and break-up of the initial sheets  and the following layer-to-layer interaction of the resulting fragments finally produced a curved and disordered flake with a number of stacked layers that, at some places, was larger than the initial figure of six. Figure 5 shows the resulting flake after being annealed by cooling it at low temperature (plots have been produced using VMD software [45] ). The top view in Figure 5(a) displays a turbostratic arrangement of the layers of material within the flake, consistent with the experimental observations. The side view of Figure 5(b) shows the curvature of the flake and the existence of stacked layers that can clearly exceed the number of six even at the edge of the flake. There are portions of the flake interior that display a larger numbers of layers. Figure 6 shows the radial distribution function (RDF) of the atoms in the flake processed at 900 K. The RDF of a single sheet of material at very low temperature is also plotted for comparison purposes. This second curve has been arbitrarily scaled to facilitate its contrast. The main peaks due to the in-layer interatomic distances in the single sheet of material are also clearly distinguishable in the sample that has been processed at high temperature. Nevertheless, the disorder in the flake is evident from the broadening of these peaks. In the figure, the peak at 3.21 A in the signal corresponding to the sample heated to 900 K has been highlighted. Even though there are nearby contributions from in-layer distances, the height of this peak relative to that of the low temperature single-sheet curve is much higher than in the neighbouring peaks. Therefore, this feature can be predominantly attributed to the interlayer stacking. The separation obtained from the calculations was in excellent agreement with the measured values at the corresponding processing temperature.

X-ray diffractometry (XRD), scanning (SEM), and transmission electron microscopy (TEM) observations
A marked corrugation of the polymeric carbon sheets, synthesized from melamine cyanurate (discussed in [19] ), was confirmed at temperatures above 600 K. Particles with tubular shape, formed by crumpling of flake particles, were observed at synthesis temperatures around 700 K (Fig. 7).
In fact, at 900 K, the crumpled particles, formed not only tubular particles but also globular ones, as tubules were further shrunk up to form globules with holes (Fig. 7).
This characteristic should be linked to H-bonding, and -as simulations showed-borders would be more instable and tend to crumple, transforming the flakes in either tubular or globular particles.
It is noteworthy to point out that this folding effect is coupled to a tighter interlayer distance upon synthesis temperature increase, indicating that -concomitant to the borders' foldingthe stacking order increases, as shown in Figure 8. This phenomenon is in agreement with the peak found at 3.21 A in RDF. Therefore, the increase in the number of stacked sheets and the larger interlayer interactions would be linked to the folding effect that produces crumpled particles. In a previous research, it was found that the layers behaved as independent ones (i.e., as a 2D material) below a synthesis temperature of 800 K (around 650 C). Above this temperature, the folding effect increased the number of stacked layers as well as the order along the c direction. In summary, the diminution of the interaction on the layer plane implicates a stronger interplanar Hbond interaction.
This corrugation effect could be conveniently used to produce globular, spherical or tubular particles. Moreover, this effect of folding seemed to affect THz absorption in a considerable way, [19] so that it could be used as a sensor of pressure in special environments. Hydrogen bond in carbon nitride (obtained from urea) was experimentally investigated by Hu et al. [46] and they concluded that melem oligomers were formed and arranged in layers. In spite of that study, we suggest that   Figure 5 and for a single sheet of material at low temperature.
the polymeric structure is compatible with the crumpled morphology. The weakening of hydrogen bonding between adjacent polymer ribbons is the main cause of the typical crumpled structure of this material, especially observed at high synthesis temperatures (ca. 900 K). This process makes the material more sensitive to the external environment through the possibility of interacting by hydrogen bonding with molecules such as CO 2 , etc.

Conclusion
A computational study of the morphological properties of a carbon nitride polymer was performed at the nanoscale using molecular dynamics. In spite of the fact that this survey was restrained by the limited size of the simulation domains, a good agreement with the experimental results was found. The synthesized samples consisted of crumpled nanosheets formed by a turbostratic arrangement of polymer layers, which is consistent with the small-size molecular dynamics simulations at the corresponding processing temperatures.
The analysis of the computer simulations results have provided a clear picture of the mechanisms that would lead to such arrangement. One determining factor is the weakness of the hydrogen bonding network that keeps together the 2D polymer and that determines many of its physical and chemical properties. This fact had been previously noted, for instance, in relation with the calculation of the electronic band-gap of this material. Due to the very small interactions between heptazine units, the calculated bandgap for the 2D system is very similar to that of the monomer. At sufficiently high temperatures, molecular agitation can partly release this relatively weak bonding network and permit the conformational change in the material. The second element is the effect of the layer-to-layer interaction, Figure 7. TEM and SEM images of polymeric carbon nitride obtained at 700 and 900 K from melamine cyanurate pyrolysis. Figure 8. Variation of the interlayer distance and stacked sheets as a function of the synthesis temperature. [19] typical of 2D carbon materials, that drives such change. Hydrogen bonding weakening between adjacent polymer ribbons is the main cause of the typical crumpled structure of this material, especially observed at synthesis temperatures above ca. 900 K. However, this property makes this material more sensitive to the external environment, and thus more suitable to be used as a catalyst and a sensor.