Ab-Initio Molecular Dynamics Simulation of Graphene Sheet

The study of graphene is important because it is a promising material for a variety of applications in the electronic industry. In the present work, the properties of а 2D periodic graphene sheet are studied with the use of ab initio molecular dynamics. DFT in the generalized gradient approximation is used in order to carry out the dynamical simulations. The PBE functional and DZVP-MOLOPT basis set are implemented in the CP2K/Quickstep package. A periodic box, consisting of 288 carbon atoms is chosen for the simulations. After geometry optimization it has dimensions 2964 x 2964 x 1500 pm and form angles of 90, 90, 60 degrees. The dynamical simulation is run for 1 ps in the NPT ensemble, at temperature T = 298.15 K. The radial distribution function shows a first peak at 142 pm, marking the bond length between carbon atoms. The density of states for the periodic systems is simulated as occupied orbitals represent the valence band and unoccupied ones the conduction band. The calculated bandgap, as expected is close to 0 eV.


Introduction
The study of graphene is important as it is a promising new material for a variety of applications in the electronic industry, as well as a highly durable or low resistant coating and medical applications [1]. Theoretical and computational approaches in studying conjugated bonds' surfaces such as graphene include mainly classical molecular dynamics and ab initio methods such as Density Functional Theory (DFT) and Post-Hartree-Fock methods. Classical molecular dynamics allow for exploring larger systems, but generally suffer in describing the electronic structure of the material as it relies on predefined force fields and potentials to represent bonding between atoms. Тhe phonon properties of graphene have been investigated using classical molecular dynamics with Tersoff-2010, LCBOP28 and AIREBO29 potentials [2]. Molecular dynamics also allows the exploration of hybrid systems where graphene interacts with polymers [3].
Electronic effects in graphene are theoretically studied using ab initio static geometry optimization methods that allow electronic structure calculations. Density Functional Theory in generalized gradient approximation (GGA) is employed in exploring doped and defected structures [4]. Further, ab initio molecular dynamics offers a serious advantage over static geometry optimization methods: it allows time dependent investigation of the system's properties. For example, the Born-Oppenheimer molecular dynamics with GGA functional was used to study the properties of a graphene sheet using 60 carbon atoms in a periodic system and the most important structural factors were revealed, even with such a relatively small system [5].
As we began a program for a thorough study of ordered and disordered carbon including simulations of the density of states, vibrational frequencies etc., we aimed at several goals: (i) to conduct ab initio simulations using a large cell, (ii) to optimize the interatomic distances in the elementary cell and (iii) to verify our results for the molecular orbitals of graphene and compare them with experimental data. Here we present some preliminary results of ab initio simulations using a larger periodic system consisting of 288 carbon atoms.

Computational details
Ab Initio Born -Oppenheimer molecular dynamic simulations were performed using the CP2K/Quickstep package [6,7]. DFT was applied within GGA, using Perdew-Burke-Ernzerhof (PBE) functional [8]. The basis set DZVP-MOLOPT-SR-GTH [9], which is optimized for calculating molecular properties in gas and condensed phase, was applied for all atoms in the studied systems. For reducing computational cost Gaussian and Plane-Wave (GPW) method was used [10]. This method uses an atom-centered Gaussiantype basis to describe the wave functions and an auxiliary plane wave basis to describe the electron density [11]. It is worth noting that only the valence electrons are explicitly treated. Their interaction with the remaining ions is described using the pseudopotentials of Goedecker-Teter-Hutter (GTH) [12,13].
The charge density cutoff of the finest grid level is equal to 400 Ry. The number of used multigrids is 5. All Born -Oppenheimer molecular dynamic simulations were carried out in the NPT ensemble with a timestep of 0.5 fs. The temperature was set to constant 298.15 K using a canonical sampling through velocity rescaling (CSVR) thermostat. Statistics were collected after thermal equilibration. The studied system consisted of a single graphene sheet represented by 288 carbon atoms.

Results and discussion
The elementary cell of graphene, used to prepare the periodic system is presented in figure 1. The triclinic cell consists of two carbon atoms and has dimensions (a = b) of 246 pm and the angle between C-C bonds is γ = 120 0 . Such a triclinic cell is chosen over orthogonal as it permits for proper quantization of the Brillouin zone, proper description of permitted k-vectors and thus offers a realistic description of the structural and electronic properties of the system [14]. The supercell used for dynamical simulations consists of a periodic system of 288 carbon atoms that occupy a hexagonal cell with dimensions A=12a and B=12b (i.e. 144 duplications of the elementary cell of graphene in the basis of the supercell) while the C-edge was fixed to 1500 pm. The volume along the C-edge above the basic graphene layer remains unoccupied. The dimension of the C-edge of the supercell is assumed to be long enough to eliminate the interaction between the graphene layers in the periodic system. The angles between A, B and C edges of the supercell are as follows: α-between A and C, β-between B and C and γ-between A and B with values of 90 0 , 90 0 and 60 0 , respectively. In order to prepare the initial system, geometry optimization is performed and the energy of the system is optimized taking into account the positions of the atoms as variables, e.g. the optimization is not performed on a structure with predefined overall dimensions but it is optimized computationally as well. After the geometry optimization, the dimensions obtained are: A= 2964 pm, B= 2964 pm and C= 1500 pm. The system should behave as a single layer. The actual system after optimization is presented in figure 2. Born -Oppenheimer molecular dynamics of the optimized system was continued for 1 ps. The NPT ensemble was employed at pressure 1 atm, in order for the size of the system to be able to change with the vibrational and dynamical changes. Thermal equilibration is observed after 200 fs. After 1 ps, the radial distribution function (RDF) of one carbon atom versus the whole system was plot, the RDF is presented in figure 3. The first peak (at 142 pm) corresponds to the 3 neighboring C atoms of the first coordination shell, e.g. it is the actual C-C interatomic distance. The second peak at 246 pm corresponds to the 6 neighboring C atoms from the second coordination shell and the third peak at 286 pm corresponds to the 3 neighboring C atoms of the third coordination shell. The experimentally measured by transmission electron microscopy average C-C bond length in single-crystalline regions of single-layered graphene is also 142 pm, with a standard deviation of 5 pm [15,16]. The theoretical data [17] suggest that graphene does not have a bandgap as the conduction and valence bands coincide at the Dirac point. The properties of graphene depend on dimensions, intrinsic defects, even on the boundaries of the specimens studied. Ab initio calculations [18] show that energy gaps arise at the edges of graphene nanoribbons (GNR). Recently published experimental results show that: i) GNR can have a bandgap of about 1.5 eV and form heterojunctions [19] and ii) confirm that the zigzag edges of graphene nanoribbon open up a bandgap of 1.9 eV at the zigzag edge states [20]. The results obtained in our study, however, do not depend on the edge states effect, as periodic boundary conditions that treat the system as an infinite sheet with no edges are employed. Additionally, the system is sufficiently large (288 atoms) to represent correctly (in sufficient amount) the delocalized π and π* orbitals around HOMO and LUMO that form the valence and conduction bands.

Conclusions
The Born -Oppenheimer molecular dynamics with the Density Functional PBE represents correctly the structural properties of a single graphene layer and especially the C-C bond lengths. The optimized bond length is 1.42 pm which exactly coincides with the experimental findings. The electronic structure and density of states were also represented correctly and a bandgap close to 0 eV was obtained. The possibility to examine complex 2D systems with GGA Density Functional Theory provides a prospective tool for exploring functionalized or multilayered graphene as well as different sp 2 /sp 3 mixed carbon phases.