Periodic DFTB for Supported Clusters: Implementation and Application on Benzene Dimers Deposited on Graphene

Abstract

The interest for properties of clusters deposited on surfaces has grown in recent years. In this framework, the Density Functional based Tight Binding (DFTB) method appears as a promising tool due to its ability to treat extended systems at the quantum level with a low computational cost. We report the implementation of periodic boundary conditions for DFTB within the deMonNano code with k-points formalism and corrections for intermolecular interactions. The quality of DFTB results is evaluated by comparison with dispersion-corrected DFT calculations. Optimized lattice properties for a graphene sheet and graphite bulk are in agreement with reference data. The deposition of both benzene monomer and dimers on graphene are investigated and the observed trends are similar at the DFT and DFTB levels. Moreover, interaction energies are of similar orders of magnitude for these two levels of calculation. This study has evidenced the high stability of a structure made of two benzene molecules deposited close to each other on the graphene sheet. This work demonstrates the ability of the new implementation to investigate surface-deposited molecular clusters properties.

1. Introduction

The modeling of functional extended surfaces has grown in past decades to investigate, for fundamental and engineering purposes, a large number of phenomena or applications such as, e.g., deposition [1], growth and migration [2], 2D assembly [3], catalysis [4], electrocatalysis [5], photocatalysis [6], molecular electronics [7], depollution [8], sensing [9], etc. Many of these studies have focused on deposited clusters, i.e., finite aggregations of basis elements (atoms or molecules) adsorbed on surfaces. Indeed, the physico-chemical properties of a cluster are distinct from the ones of both the single entities and the infinite cluster (bulk) and strongly depend on the size and structure of the cluster [10]. Understanding and controlling the structure of deposited clusters could thus allow the precise tuning of their properties.

The theoretical study of clusters deposited on extended surfaces is very challenging due to the size of the space of structural and electronic configurations to be explored and the high level of computational methods that has to be implemented. Indeed, the a priori unknown nature of the interactions between the cluster building blocks and between the cluster and the surface (with potential charge transfers at the cluster–surface interface) prevents the use of empirical force fields in favor of methods in which the electronic structure is explicitly considered. The very high computational cost of ab initio and Density Functional Theory (DFT) methods prohibit their use for the study of such systems, particularly when dealing with global optimisation or finite-temperature molecular dynamics, and one strategy consists in implementing, in a periodic formulation, approximate quantum mechanical methods. Among them, DFTB is an approximated DFT scheme with a much lower computational cost enabled by the use of parameterized integrals in a minimal valence basis set [11,12,13,14]. There have been several implementations of DFTB within periodic conditions in various codes [15,16,17,18,19,20,21], allowing the computation of structural, mechanical and electronic properties. In the present work, we report a new implementation of the DFTB scheme for periodic systems within the deMonNano code [22] combined with corrections to describe long range intermolecular interactions.

After assessing the performance of this implementation on graphene and graphite, we apply this method to characterize model systems consisting in benzene and benzene dimers deposited on a graphene surface. Such a system is relevant in an astrophysical context as it can be seen as a simple model of Polycyclic Aromatic Hydrocarbons (PAH) clusters adsorbed on large carbonaceous grains or on very large PAHs. Indeed, despite the fact that PAHs are expected to be ubiquituous in the interstellar medium [23] and their clusters have been proposed to play a significant role in interstellar physics and chemistry [24], the structural and energetic property changes induced by their deposition on a surface remain, to a large extent, unknown. The second motivation for selecting this benchmark system is that a reasonable description of the benzene dimers potential energy surface is challenging even with ab initio schemes [25], making it a system of choice to address the quality of our approach. This is due to the fine equilibrium between Pauli repulsion, dispersion and coulomb interaction, which drives the competition between parallel and T-shaped structures. In the past, we have shown that the combination of DFTB with empirical dispersion and atomic charges corrections allowed for a proper description of such systems [26].

In this paper, the periodic formulation of DFTB that has been implemented in deMonNano is presented in Section 2, with a special focus on the originality of the present scheme with respect to other periodic implementations, i.e., its combination with the WMull charge-correction approach. Computational details are provided in Section 3, and the benchmarks on graphene and graphite are presented in Section 4. The applications to benzene monomer and dimers deposited on graphene are discussed in Section 5. Finally, a conclusion is provided in Section 6.

2. Methods

2.1. DFTB

The Density Functional-based Tight-Binding method (DFTB) can be derived from DFT from several approximations [11,13,14,16,27].

The first one relies on the expression of molecular orbitals (MOs) ${\varphi }_i\left(\mathbf{r}\right)$ as linear combinations of atomic orbital (LCAO)-type basis sets using minimal valence bases ${\chi }_{\mu }$:

$${\varphi }_i\left(\mathbf{r}\right)=\sum _{\mu }{c}_{i\mu }{\chi }_{\mu }\left(\mathbf{r}\right)$$

(1)

where $c_{i\mu }$ is the coefficient of the molecular orbital i on the atomic orbital ${\chi }_{\mu }$.

A Taylor expansion of the DFT energy is performed as a function of the electronic density, and the real density $\rho$ of the system minimizing the Kohn–Sham energy is searched as a perturbation with respect to reference density ${\rho }_0$ ($\rho ={\rho }_0+\delta \rho$).

$$\begin{array}{cc}\hfill E\left[\rho \left(\mathbf{r}\right)\right]=E\left[{\rho }_0\left(\mathbf{r}\right)\right]+& \int {\left.\frac{\delta E\left[\rho \right(\mathbf{r}\left)\right]}{\delta \rho \left(\mathbf{r}\right)}\right|}_{{\rho }_0}\delta \rho \left(\mathbf{r}\right)+\frac{1}{2}\int \int {\left.\frac{{\delta }^{2}E\left[\rho \left(\mathbf{r}\right)\right]}{\delta \rho \left(\mathbf{r}\right)\delta \rho \left({\mathbf{r}}^{\prime }\right)}\right|}_{{\rho }_0}\delta \rho \left(\mathbf{r}\right)\delta \rho \left({\mathbf{r}}^{\prime }\right)+\hfill \\ \hfill & \cdots +\frac{1}{p!}\int \int \cdots \int {\left.\frac{{\delta }^{p}E\left[\rho \left(\mathbf{r}\right)\right]}{\delta \rho \left(\mathbf{r}\right)\delta \rho \left({\mathbf{r}}^{\prime }\right)\cdots \delta \rho \left({\mathbf{r}}^{\left(\mathbf{p}\right)}\right)}\right|}_{{\rho }_0}\delta \rho \left(\mathbf{r}\right)\delta \rho \left({\mathbf{r}}^{\prime }\right)\cdots .\delta \rho \left({\mathbf{r}}^{\left(\mathbf{p}\right)}\right)\hfill \end{array}$$

(2)

In the original version, also known as the non self-consistent DFTB (sometimes referred to as zeroth-order DFTB or simply DFTB [11,12]), only the zeroth and first-order terms of the Taylor expansion are retained. In the DFTB2 scheme [27], also known as self-consistent charge (SCC) DFTB and in the DFTB3 scheme [28], the second-order and third-order terms are also taken into account, respectively.

At the DFTB0 level, the electronic energy reads as follows:

$$E^{DFTB0}=\sum _{\alpha <\beta }{E}_{rep}\left({\mathbf{r}}_{\alpha \beta }\right)+\sum _{i\mu \nu }{n}_i{c}_{i\mu }{c}_{i\nu }{H}_{\mu \nu }^0$$

(3)

with ${\mathbf{r}}_{\alpha \beta }={\mathbf{r}}_{\beta }-{\mathbf{r}}_{\alpha }$, where ${\mathbf{r}}_{\alpha }$ and ${\mathbf{r}}_{\beta }$ are the positions of atoms $\alpha$ and $\beta$, respectively, $E_{rep}\left({\mathbf{r}}_{\alpha \beta }\right)$ is a repulsive contribution between atoms $\alpha$ and $\beta$, $n_i$ is the occupation of the orbital i and $H_{\mu \nu }^0$ is the matrix element associated with the Kohn–Sham operator at the reference density expressed in the atomic basis. These matrix elements, as well as those of the atomic basis overlap matrix S, can be parameterized as only one-body or two-body terms. This is allowed by the definition of the reference density as a superposition of atomic densities ${\rho }_0={\rho }_0^{\alpha }+{\rho }_0^{\beta }+{\rho }_0^{\gamma }+\dots$, and the reduction of integrals to one-center or two-center terms:

• $H_{\mu ,\nu \in \alpha }^0\left({\rho }_0\right)\approx H_{\mu \nu }^0\left({\rho }_0^{\alpha }\right)\approx {\delta }_{\mu \nu }{ϵ}_{\mu \alpha }$; the atomic orbital energies of the isolated atom $\alpha$;
• $H_{\mu \in \alpha ,\nu \in \beta }^0\left({\rho }_0\right)\approx H_{\mu \nu }^0({\rho }_0^{\alpha }+{\rho }_0^{\beta })$, which only depends on the distance between the two corresponding atomic centers: $H_{\mu \in \alpha ,\nu \in \beta }^0({\mathbf{r}}_{\alpha }-{\mathbf{r}}_{\beta })$.

Focusing from now on the SCC-DFTB level [27], the energy expression is:

$$E^{SCC-DFTB}=E^{DFTB0}+\frac{1}{2}\sum _{\alpha \beta }{\gamma }_{\alpha \beta }{q}_{\alpha }{q}_{\beta }$$

(4)

The last term corresponds to the second-order contribution and depends on the electronic density fluctuation $\delta \rho$ represented by atomic charges $q_{\alpha }$. ${\gamma }_{\alpha \beta }$ is a matrix for which diagonal terms are equal to the atomic Hubbard parameters and off-diagonal terms contain the $1/R$ coulomb interaction between atomic charges and an exchange-correlation energy $E_{xc}$ contribution:

$${\gamma }_{\alpha \beta }=\int \int \left(\frac{1}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}+{\left.\frac{{\delta }^2{E}_{xc}}{\delta \rho \left(\mathbf{r}\right)\delta \rho \left({\mathbf{r}}^{\prime }\right)}\right|}_{{\rho }_0}\right)F_0^{\alpha }(\mathbf{r}-{\mathbf{r}}_{\alpha })F_0^{\beta }({\mathbf{r}}^{\prime }-{\mathbf{r}}_{\beta })d\mathbf{r}d{\mathbf{r}}^{\prime }$$

(5)

where $F_0^{\alpha }$ is the normalized spatial extension for the excess/default of electrons around atom $\alpha$ with respect to the neutral atom, and it is assumed here to have no angular dependence.

Since the second-order term contains atomic charges, this introduces a term depending on the charges of $H^1\left(q\right)$ into the TB operator:

$$(H^0+H^1\left(q\right))C_i={ϵ}_{i}S{C}_i$$

(6)

with the following:

$$H_{\mu \nu }^1=\frac{1}{2}{S}_{\mu \nu }\sum _{\xi }{q}_{\xi }({\gamma }_{\alpha \xi }+{\gamma }_{\beta \xi })$$

where $\mu$ and $\nu$ belong to atoms $\alpha$ and $\beta$, respectively, and $C_i$ is the column vector containing the coefficients for the ith MO. As the charges depend on MO coefficients $c_{i\mu }$, the new secular equation must be solved self-consistently with respect to atomic charges, at the origin of the method’s name self-consistent charge (SCC-)DFTB.

In the standard SCC-DFTB version [27], the atomic charges are computed from the density matrix P and the atomic basis overlap S matrix within the Mulliken approximation:

$$q_{\alpha }=\sum _{\mu \in \alpha }\sum _{\beta }\sum _{\nu \in \beta }{P}_{\mu \nu }{S}_{\mu \nu }$$

(7)

making use of the density matrix:

$$P_{\mu \nu }=\sum _{\mu \nu }{n}_i{c}_{i\mu }{c}_{i\nu }$$

In previous works, we have shown that atomic charges can be improved by taking into account the bond polarisation, adapting the Charge Model class IV scheme for DFTB [26,29]. This approach, which requires the calculation of Mayer’s bond order, is computationally expensive and hardly transferable to a periodic implementation. To circumvent this bottleneck, we have recently introduced a simpler scheme, named in the following WMull for Weighted Mulliken charges [30], to correct atomic charges with the following expression.

$$q_{\alpha }=\sum _{\mu \in \alpha }\sum _{\beta }\sum _{\nu \in \beta }{P}_{\mu \nu }{S}_{\mu \nu }(1+t_{\alpha \beta })$$

(8)

$t_{\alpha \beta }=-t_{\beta \alpha }$ is an empirical parameter accounting for a non-symmetric repartition of the electrons between different atomic types, where the Mulliken symmetric repartition is recovered for $t_{\alpha \beta }=0$. The second order contribution to the Kohn–Sham operator matrix is modified as follows.

$$H_{\mu \nu }^1=\frac{1}{2}{S}_{\mu \nu }\sum _{\xi }{q}_{\xi }({\gamma }_{\alpha \xi }(1+t_{\alpha \beta })+{\gamma }_{\beta \xi }(1-t_{\alpha \beta }))$$

(9)

We have shown that this simple scheme provides similar results to those obtained with the Charge Model approach to model clusters of PAHs [31] and water [32].

2.2. DFTB for Periodic Systems

The former implementation of periodic DFTB within deMonNano was restricted to $\mathrm{\Gamma }$-point approximation only [22]. In the present implementation, after defining a set of k-points in the reciprocal space, the electronic problem is searched self consistently. A step of the self-consistent scheme consists in solving separately the secular equations for each k-point to obtain molecular orbitals ${\varphi }_i^{\mathbf{k}}$. Molecular orbitals obtained for all k-point are then used to build the total electronic density, the latter being used as an input for the next self-consistent step.

For a given k-point, the molecular orbitals ${\varphi }_i^{\mathbf{k}}$ must fulfill Bloch’s theorem:

$${\widehat{T}}_{\mathbf{R}}{\varphi }_i^{\mathbf{k}}=e^{i\mathbf{kR}}{\varphi }_i^{\mathbf{k}}\left(\mathbf{r}\right)$$

(10)

where ${\widehat{T}}_{\mathbf{R}}$ is the operator associated to a translation of $\mathbf{R}$ and where $\mathbf{R}$ is a vector connecting two unit cells. This is achieved by expanding MOs on a basis of Bloch functions built from real space atomic orbitals:

$${\varphi }_i^{\mathbf{k}}\left(\mathbf{r}\right)=\sum _{\mathbf{k}}{c}_{i\mu }^{\mathbf{k}}{\chi }_{\mu }^{\mathbf{k}}\left(\mathbf{r}\right)$$

$${\chi }_{\mu }^{\mathbf{k}}\left(\mathbf{r}\right)=\frac{1}{\sqrt{N}}\sum _N{e}^{i{\mathbf{kR}}_N}{\chi }_{\mu }\left(\mathbf{r}-{\mathbf{R}}_N\right)$$

(11)

where the infinite sum relies on all N possible translations from the main unit cell to any of the other ones through translation vectors ${\mathbf{R}}_N$. The overlap and Kohn–Sham operator matrices expressed in this basis can be written from their real-space equivalent by making use of the following transformation:

$$H_{\mu \nu }^{\mathbf{k}}=\sum _N{e}^{i{\mathbf{kR}}_N}{H}_{\mu \nu }({\mathbf{r}}_{\alpha }-{\mathbf{r}}_{\beta }-{\mathbf{R}}_N)=\sum _N{e}^{i{\mathbf{kR}}_N}(H_{\mu \nu }^0({\mathbf{r}}_{\alpha }-{\mathbf{r}}_{\beta }-{\mathbf{R}}_N)+H_{\mu \nu }^1({\mathbf{r}}_{\alpha }-{\mathbf{r}}_{\beta }-{\mathbf{R}}_N))$$

and

$$S_{\mu \nu }^{\mathbf{k}}=\sum _N{e}^{i{\mathbf{kR}}_N}{S}_{\mu \nu }({\mathbf{r}}_{\alpha }-{\mathbf{r}}_{\beta }-{\mathbf{R}}_N)$$

(12)

where $\mu$ and $\nu$ belong to atoms $\alpha$ and $\beta$, respectively.

In the previous expressions, the matrix elements of $H^0$ and S are easily obtained from the DFTB Slater–Koster tables [33] and rapidly vanish for large values of ${\mathbf{R}}_N$. The first order contribution to $H_{\mu \nu }$ is also short range with respect to $({\mathbf{r}}_{\alpha }-{\mathbf{r}}_{\beta }-{\mathbf{R}}_N)$ but contains an infinite long range coulomb sum.

$$H_{\mu \nu }^{1,\mathbf{k}}({\mathbf{r}}_{\alpha }-{\mathbf{r}}_{\beta }-{\mathbf{R}}_N)=\frac{1}{2}{S}_{\mu \nu }({\mathbf{r}}_{\alpha }-{\mathbf{r}}_{\beta }-{\mathbf{R}}_N)\sum _{\xi }\sum _N{q}_{\xi }({\gamma }_{\alpha \xi }({\mathbf{r}}_{\alpha }-{\mathbf{r}}_{\xi }-{\mathbf{R}}_N)+{\gamma }_{\beta \xi }({\mathbf{r}}_{\beta }-{\mathbf{r}}_{\xi }-{\mathbf{R}}_N))$$

(13)

In practice, this infinite sum is replaced by an Ewald summation. The secular equation is solved for each k-point.

$$H^{\mathbf{k}}{C}_i^{\mathbf{k}}={ϵ}_i^{\mathbf{k}}{S}^{\mathbf{k}}{C}_i^{\mathbf{k}}$$

(14)

The eigenvalues ${ϵ}_i^{\mathbf{k}}$ resulting from all the k-point secular equations are then sorted in ascending order to drive the determination of the orbital occupation numbers $n_i^{\mathbf{k}}$ following either a canonical occupation or a Fermi distribution.

The density matrix can, therefore, be computed for each k-point.

$$P_{\mu \nu }\left(\mathbf{k}\right)=\sum _i{n}_i^{\mathbf{k}}{c}_{i\mu }^{\mathbf{k}\ast }{c}_{i\nu }^{\mathbf{k}}$$

(15)

We follow the approach of reference [15], which consists in building the real space density matrix and computing atomic charges in the real space. The real space density matrix is obtained by summing over k-points:

$$P_{\mu \nu }\left({\mathbf{R}}_N\right)=\sum _{\mathbf{k}}{P}_{\mu \nu }\left(\mathbf{k}\right)e^{-i{\mathbf{kR}}_N}$$

(16)

and the Mulliken charges are then computed as follows:

$$q_{\alpha }=\sum _{R_N}{P}_{\mu \nu }\left({\mathbf{R}}_N\right)S_{\mu \nu }({\mathbf{r}}_{\alpha }-{\mathbf{r}}_{\beta }-{\mathbf{R}}_N)$$

(17)

and they are used as inputs for the next SCC cycle.

The simple WMull correction to Mulliken charges can be generalized to the periodic equations replacing Equation (17) by the following:

$$q_{\alpha }=\sum _{R_N}{P}_{\mu \nu }\left({\mathbf{R}}_N\right)S_{\mu \nu }({\mathbf{r}}_{\alpha }-{\mathbf{r}}_{\beta }-{\mathbf{R}}_N)(1+t_{\alpha \beta })$$

(18)

and Equation (13) by the following.

$$\begin{array}{ccc}\hfill H_{\mu \nu }^{1,\mathbf{k}}({\mathbf{r}}_{\alpha }-{\mathbf{r}}_{\beta }-{\mathbf{R}}_N)=\frac{1}{2}{S}_{\mu \nu }({\mathbf{r}}_{\alpha }-{\mathbf{r}}_{\beta }-{\mathbf{R}}_N)\sum _{\xi }\sum _N& q_{\xi }& ({\gamma }_{\alpha \xi }({\mathbf{r}}_{\alpha }-{\mathbf{r}}_{\xi }-{\mathbf{R}}_N)(1+t_{\alpha \beta })\hfill \\ & +& {\gamma }_{\beta \xi }({\mathbf{r}}_{\beta }-{\mathbf{r}}_{\xi }-{\mathbf{R}}_N)(1-t_{\alpha \beta }))\hfill \end{array}$$

(19)

3. Computational Details

3.1. DFTB Calculations

Different DFTB parameters are available in the literature (cf. website www.dftb.org) depending on the choices made during the parameterization procedure such as the DFT functional, the basis sets type (Gaussian and Lorentzian) used to generate the atomic orbitals, the confinement imposed on these orbitals, the reference data used to compute the repulsive contribution E${}_{rep}$ and, for the second and third order DFTB, the values of the atomic Hubbard parameters and their derivatives. In this work, we are working with the BIO DFTB set of parameters [27] provided within the deMonNano code.

Dispersion interaction corrections can be introduced in the DFTB Hamiltonian using an empirical atom-atom formula. Two types of corrections are available in the deMonNano code and will be tested in the next section. The first one (hereafter labelled D1) is a Lennard–Jones type potential with short range corrections introduced by Zhechkov et al. [34]. The second one (hereafter labelled D2 [26]) is given by the following expression:

$$E_{\mathrm{disp}}=-\sum _N\sum _{\alpha ,\beta }{f}_{\mathrm{damp}}\left(\right|{\mathbf{r}}_{\alpha }-{\mathbf{r}}_{\beta }-{\mathbf{R}}_N\left|\right)\frac{C_{\alpha \beta }^6}{|{\mathbf{r}}_{\alpha }-{\mathbf{r}}_{\beta }-{\mathbf{R}}_N{|}^6}$$

(20)

where $f_{\mathrm{damp}}$ is a damping function screening the short range contribution and $C_{\alpha \beta }^6$ is an empirical parameter (see [26] for details). In both cases, only the van der Waals contributions larger than ${10}^{-5}$ Hartree are taken into account in order to limit the number N of boxes involved in the sum.

When calculations are performed with the WMull scheme, a value of $t_{CH}=0.245$ has been determined to provide the atomic charges for the benzene molecule in agreement with reference calculations (see tables and discussion in reference [26]).

Regarding convergency criterions, we have used a tolerance of ${10}^{-8}$ for the atomic charges during the SCC process and 5 $\times {10}^{-6}$ Hartree/Bohr for the largest gradient for local optimizations.

3.2. DFT Calculations

Dispersion-corrected DFT calculations were performed under periodic boundary conditions using the Vienna ab initio simulation package (VASP [35,36,37]) together with PAW pseudopotentials [38,39] and the PBE-D3 semiempirical dispersion-corrected functional in its zero-damping formalism [40]. This functional has been chosen as it has been reported as a relevant choice for studies involving graphene [41]. A conjugate-gradient algorithm was used to relax the ions and the convergence criterion was set up so that the maximum atomic force was less than 0.01 eV Å${}^{-1}$, and all atoms were allowed to relax unconstrained. To avoid interactions between the benzene monomers/dimers and their periodic images, a cubic box measuring 50 Å on a side was used for isolated systems. For supported ones, the calculations were performed on a 29.92 Å × 34.55 Å graphene surface (these values having been calculated on the basis of the graphene equilibrium lattice parameter reported in Section 4) containing 392 carbon atoms placed in a 50 Å high simulation box to avoid any interaction between the adsorbed molecule and the underside of the graphene sheet of the upper periodic box. Since the size of the supercell was large enough, the Brillouin zone sampling in reciprocal space restricted to the $\mathrm{\Gamma }$-point was sufficient to ensure good convergence of the total energy, except for the calculations aimed at determining the equilibrium parameters of the graphite bulk, which required a 1 × 1 × 5 k-points grid. A plane-wave kinetic energy cutoff of 450 eV was employed. For dealing with the partial occupancies around the Fermi level, Methfessel–Paxton smearing was used with $\sigma$ = 0.2 eV [42].

4. Benchmark Calculations: Graphene and Graphite

The goal of this section is twofold. First, k-points grid size tests are performed on two systems of interest, namely the 392-atom graphene sheet that will be further used in Section 5 and 2 × 392-atom graphite built from periodic repetition of two graphene layers. Second, we aim at evaluating the quality of the DFTB potential to reproduce the main features of graphene and graphite. In order to model graphene, we have first optimized the DFTB lattice parameter, working with a periodic box containing 392 atoms (∼30 Å × 35 Å × 50 Å). For such a large simulation box, the $\mathrm{\Gamma }$-point approximation remains valid as the energy varies by less than 3.4 $\times {10}^{-4}$ eV/atom (0.0008%) when increasing from one to three k-points in the x and y directions and by less than 1.3 $\times {10}^{-4}$ eV/atom (0.0003%) when changing from three to five k-points. The equilibrium C-C bond lengths determined with one or three k-points in x and y directions are the same at the precision of ${10}^{-3}$ Å. Values of 1.430 ± 0.001 Å and 1.426 ± 0.001 Å were obtained with DFTB-D1 and DFTB-D2 methods, respectively (see

Table 1. Graphene and Graphite structural data (in Å). * In these computational studies, d${}_{C-C}$ distance was fixed to the one determined experimentally. ** Differences between the two DFTB-D1 calculations are detailed in the text.

Methods	d${}_{\mathit{C}-\mathit{C}}^{\mathit{Graphene}}$	d${}_{\mathit{C}-\mathit{C}}^{\mathit{Graphite}}$	d${}_{\mathit{interlayer}}^{\mathit{Graphite}}$
DFTB-D1 **	1.430 ± 0.001	1.429 ± 0.001	3.383 ± 0.001
DFTB-D2	1.426 ± 0.001	1.425 ± 0.001	3.131 ± 0.001
DFTB-D1 **	1.421 [34]	1.421 [34]	3.38 [34]
DFT-D3	1.425 ± 0.001	1.424 ± 0.001	3.488 ± 0.001
PW91-OBS	1.4226 [43]
RPA		1.42 * [44]	3.34 * [44]
QMC		1.42 * [45]	3.426 * [45]
Expt	1.42 [46]	1.422 [47,48]	3.356 [47,48]

). These values are slightly larger than the value of 1.421 Å previously reported by Zhechkov et al. using the $\mathrm{\Gamma }$-point approximation and a smaller unit cell [34]. The DFTB-D2 values produces the best agreement with the C-C bond length obtained at the dispersion-corrected DFT level (1.425 ± 0.001 Å with DFT-D3 and 1.4226 Å with PW91-OBS [43]) as well as with the experimental values (1.42 Å).

The graphite bulk has been modeled by including two layers of the previously defined graphene sheet in the periodic box. In order to determine the appropriate number of k-points in the z direction (perpendicular to the graphene planes), we have performed DFTB single point energy calculations for an interlayer distance of 3.5 Å, chosen because it corresponds to the DFT-D3 one (3.488 Å, see Table 1), with one k-point in the x and y directions. The total energy varies by 2 $\times {10}^{-4}$ eV/atom (0.0004%) when the number of k-points increased from one to three in the z direction, and then it remains constant for calculations performed with five, seven and nine k-points in the z-direction. We have determined the equilibrium parameters of the graphite bulk with three k-points in the z directions (Table 1). Using either one or three k-points in the x and y directions led to the sames results at the target precision of 0.001 Å. DFTB-D1 and DFTB-D2 C-C bond lengths are reduced by 0.001 Å with respect to their values in the graphene sheet, which is a trend also observed at the DFT-D3 level. The DFTB-D1 interlayer equilibrium distance (3.383 ± 0.001 Å) is in agreement with both the value of reference [34] with a four layers model in the $\mathrm{\Gamma }$-point approximation (3.38 Å) and the experimental data (3.356 Å). The interlayer distance is reduced to 3.131 ± 0.001 Å at the DFTB-D2 level. With respect to theoretical references (DFT-D3, RPA and QMC) and experimental values, we can conclude that the DFTB-D1 method produces better quality results for graphite interlayer distances, while the DFTB-D2 method prevails for the C-C bond length.

5. Graphene Supported Benzene and Benzene Dimers

5.1. Benzene Supported on Graphene

We performed local structural optimization for systems consisting of an isolated benzene molecule deposited on top of a graphene monolayer. On the basis of the results obtained in Section 4, the calculations have been performed in the $\mathrm{\Gamma }$-point approximation, with the initial structures corresponding to a benzene molecule deposited in the proper orientation on the optimized graphene layer. Four different configurations have been probed and labeled $a_1,a_2,a_3$ and $a_4$, which can be visualized in Figure 1. The three first ones correspond to the structures labeled $a_1,a_2$ and $a_3$ in reference [49]: hollow, bridge and top in reference [50] and AA, SP and AB in reference [51]. The last structure $a_4$ was named top-rot in reference [50] and was also studied in reference [52].

The binding energies of benzene on graphene, reported on

Table 2. Binding energies of benzene on graphene in eV. The experimental binding energy of a benzene molecule on a graphite surface is −0.50 ± 0.08 eV [53].

Method	a1	a2	a3	a4
DFTB-D1	−0.639	−0.652	−0.654	−0.651
DFTB-D2	−0.439	−0.448	−0.447	−0.451
DFT-D3	−0.428	−0.450	−0.453	−0.450
LDA [49]	−0.16	−0.23	−0.24
$\mathbf{\omega }$B97X-D [51]			−0.47
optB86b-vdw [52]				−0.5
vdW-DF1 [50]			−0.49
vdW-DF2 [50]			−0.43
Expt. Saturated Adsorption Enthalpy [52]				−0.5

, were computed with the following formula:

$$E_{bind}^{graphene@bz}=E_{tot}^{graphene@bz}-E_{tot}^{graphene}-E_{tot}^{bz}$$

(21)

where all energies correspond to relaxed structures. DFTB-D1 and DFTB-D2 results agree on the main trends which are also present at the DFT-D3 level: three almost degenerated structures, namely $a_2,a_3$ and $a_4$, and the $a_1$ structure being less stable by about 0.012(DFTB-D2)/0.015(DFTB-D1)/0.025(DFT-D3) eV (see Table 2). The absolute binding energies provided by the DFTB-D2 scheme are in very good agreement with DFT results (apart from LDA) and experimental measurements. The DFTB-D1 scheme provides poorer results, with an overestimation of the binding energies of about 35% (∼0.2 eV).

Regarding z-separation between the benzene monomer and the graphene sheet (see

Table 3. Z-separation of benzene on graphene (in Å).

Method	a1	a2	a3	a4
DFTB-D1	3.15	3.13	3.13	3.14
DFTB-D2	3.08	3.07	3.08	3.05
DFT-D3	3.47	3.38	3.36	3.42
$\mathbf{\omega }$B97X-D [51]	3.36	3.30	3.35
vdW-DF1 [50]			3.6
vdW-DF2 [50]			3.5
MP2 [54]				3.33

), DFTB-D1 and DFTB-D2 benzene-graphene distances are shorter than DFT ones by ∼0.35 Å. However, it should be noted that z-separations calculated with dispersion-corrected DFT functionals also significantly differ from each other by up to 0.25 Å and experimental data are missing in the literature to our knowledge. Concerning the $a_4$ structure, z-separation can be compared with the MP2 one reported in reference [54], which confirms the slight underestimation at the DFTB-D1 and DFTB-D2 levels.

The binding energy and the z-separation of a graphene sheet with a deposited benzene monomer can be compared with the ones of larger aromatic molecules. DFTB-D2 and DFT-D3 computed benzene binding energies are in the range 0.43–0.45 eV, which is consistent with the almost three-times larger binding energies of molecules containing three aromatic rings such as anthraquinone (1.30 eV [55]) and phenanthraquinone (1.31 eV [56]). For these systems, the reported z-separations are 3.19 Å [55] for anthraquinone and 3.14 Å [56]) for phenanthraquinone, and these values last between the DFTB-D2 and DFT-D3 values obtained for the graphene-supported benzene molecule. One can also mention the slightly larger distance reported for indole (tryptophan model, 3.50 Å [54]).

5.2. Benzene Dimers in Vacuum

Reproducing the benzene dimer potential energy surface is a challenging task for DFT schemes, due to the fine competition between the various contributions to the total energy. This is even more true for approximated schemes such as the DFTB method. Briefly, three characteristic structural families can be identified, namely sandwich (S), parallel-displaced (PD) and T-shaped (T), each one presenting several minima. In this work, the sandwich-eclipsed (SE) structure has been selected to represent the sandwich family. The PD family is represented by the isomer shown in Figure 2 as it was previously reported to be the most stable of this family at the DFTB level. Two additional structures were considered to account for the T-shaped family, namely T and Csoa isomers (corresponding to T4 and Csoa in reference [26]), which only differ by a slight displacement of the top benzene from a symetric position toward a position over a carbon atom. These two structures were previously reported to be degenerated as their energies differed by less than 10${}^{-3}$ eV at the DFTB level, and the present DFT-D3 calculations show a difference of 8 $\times {10}^{-3}$ eV in favor of the Csoa isomer.

The binding energy of benzene dimer were computed with the following formula:

$$E_{bind}^{b{z}_2}=E_{tot}^{b{z}_2}-2\times E_{tot}^{bz}$$

(22)

where all energies correspond to relaxed structures. These energies are gathered in

Table 4. Binding energies for benzene dimer bz${}_2$ in eV.

Method	T-Shaped		PD	SE
DFTB-D1	Csoa	−0.126	−0.194	−0.192
DFTB-D1-WMull	Csoa	−0.140	−0.162	−0.148
DFTB-D2	T	−0.099	−0.135	−0.132
DFTB-D2-WMull	T	−0.113	−0.104	−0.086
DFT-D3	Csoa	−0.146	−0.152	−0.106
	T	−0.138
CCSD(T) [25]	Csoa	−0.12	−0.12	−0.07
SAPT [25]	Csoa	−0.12	−0.12	−0.08

for the different structures optimized with both DFTB and DFT-D3 methods, together with the ones reported in the literature for ab initio reference calculations (CCSD(T) and SAPT). In the case of the T-shaped family, only one of the two studied isomers could be located on DFTB potential energy surfaces, namely Csoa with the DFTB-D1 method and T with the DFTB-D2 one. DFT-D3 and ab initio reference calculations agree on the fact that T-shaped and PD structures are close in energy and by far more stable than the SE structure. This ordering is not reproduced at the DFTB-D1 nor at the DFTB-D2 levels, i.e., for the two dispersion corrections investigated in the absence of atomic charge corrections, because the SE structure is found to be almost degenerated with the PD structure and the T-shaped isomer is found to be the less stable one in both cases. Introducing the WMull charge correction detailed in Section 2 with the D1 dispersion (DFTB-D1-WMull) renders the PD structure the most stable, but the T-shaped structure remains the least stable. Finally, the DFTB-D2-WMull method provides a correct picture, with T-shaped and PD isomers being close in energy and more stable than the SE isomer. In addition, the binding energies are of the same order as those of the reference calculations.

5.3. Benzene Dimers Supported on Graphene

It appeared from the previous sections that the DFTB-D2-WMull method is the best DFTB choice for modeling both an isolated benzene molecule deposited on a graphene layer and a benzene dimer in vacuum. This level of theory has been chosen to investigate the deposition of a benzene dimer on a graphene sheet. Among the possible adsorption modes of a benzene molecule on a graphene monolayer, we selected the $a_4$ one as it was found to be the most stable one at this level of theory (see Section 5.1) and added a second benzene unit to form T, PD or SE configurations. The optimized structures, obtained in the $\mathrm{\Gamma }$-point approximation, are shown in Figure 2.

It can be observed that T and SE structures were preserved during optimization. On the opposite side, the deposited PD dimer was stable at the DFT-D3 level only and resulted in a dissociated configuration at the DFTB-D2-WMull level in which the two benzene units are close to each other, both exhibiting a $a_4$ adsorption configuration on the graphene sheet. This configuration is hereafter named $a_4-cd$ (close deposition on $a_4$ adsorption sites). In this latter configuration, the hydrogen atom of each benzene molecule is pointing in between two hydrogen atoms of the other benzene unit. Such a configuration limits the coulomb repulsion between the positively charged hydrogen atoms while preserving some attractive dispersion interactions.

Z-separations between the graphene sheet and the benzene units parallel to it are reported in

Table 5. Binding energy of benzene dimers bz${}_2$ on graphene in eV and z-separation between the graphene sheet and the benzene units parallel to it in Å. $*$ No value is reported at the DFTB-D2-WMull level for PD as the optimization resulted in the $a_4-cd$ structure.

${\mathit{E}}_{\mathit{tot}}^{\mathit{ref}}=$	DFTB-D2-WMull				DFT-D3
	T	PD *	${\mathit{a}}_4-\mathit{c}\mathit{d}$	SE	T	PD	${\mathit{a}}_4-\mathit{c}\mathit{d}$	SE
$E_{tot}^{graphene}+2\times E_{tot}^{bz}$	−0.586	↳	−0.929	−0.573	−0.605	−0.624	−0.938	−0.583
$E_{tot}^{graphene}+E_{tot}^{b{z}_2}$	−0.473	↳	-	−0.487	−0.467	−0.473	-	−0.478
$E_{tot}^{graphene@bz}+E_{tot}^{bz}$	−0.135	↳	−0.478	−0.122	−0.156	−0.175	−0.489	−0.134
$2\times E_{tot}^{graphene@bz}-E_{tot}^{graphene}$	0.316	↳	−0.027	0.329	0.294	0.275	−0.039	0.316
Z-separation	3.03	↳	3.00	3.02/6.53	3.38	3.39/6.91	3.40/3.40	3.39/7.33

. At both the DFT and DFTB levels, the distance between the graphene sheet and the closest benzene unit is very slightly reduced (less than 0.05 Å with DFTB-D2-WMull and 0.03 Å with DFT) with respect to the z-separation values obtained for the deposited benzene monomer (see Table 3). In the SE structure, the presence of the graphene surface results in a slight reduction in intermolecular distance from 3.55 Å to 3.52 Å at the DFTB level, whereas it slightly increases from 3.88 Å to 3.94 Å at the DFT level.

The binding energies associated with the optimized structures are computed with the following formula:

$$E_{bind}^{graphene@b{z}_2}=E_{tot}^{graphene@b{z}_2}-E_{tot}^{ref}$$

(23)

and reported in Table 5, making use of various choices for $E_{tot}^{ref}$, which is the reference potential energy. In the first line (graphene + 2bz), $E_{tot}^{ref}$ is the energy of an optimized graphene monolayer plus that of two isolated benzene molecules. It appears that, at DFT and DFTB levels, the most stable configuration relies on the dissociation of the benzene dimer to form the $a_4-cd$ structure. The energetic difference between the T-shaped structure and the less stable SE dimer is twice smaller (0.013 eV vs. 0.027 eV at the DFTB level and 0.022 eV vs. 0.052 eV at the DFT level) when the dimer is deposited with respect to the gas phase condition. This is probably due to the interaction between the graphene surface and the benzene molecule that is furthest from the surface, which is favored in the sandwich configuration. In the second line, $E_{tot}^{ref}$ is the total energy of an isolated graphene sheet plus that of the optimized dimer in its T, PD or SE form, respectively. The gained energies for the non-dissociating dimers (T and SE) are similar (∼0.46–0.49 eV for DFT and DFTB values), which is only very slightly above the binding energy of a single benzene with graphene (0.45 eV for DFT and DFTB values). In the third line, $E_{tot}^{ref}$ is the energy of a benzene deposited on a graphene sheet plus that of an isolated benzene. It differs from the isolated dimers by 0.022 eV for the T-shaped structure and 0.036 eV for the SE structure at the DFTB level and 0.018 eV and 0.028 eV at the DFT level. The larger value obtained for the SE dimer can be, again, related to the expected larger interaction energy between the graphene sheet and the second further benzene unit in the SE configuration. The last line compares the binding energies with the one of a system where two benzene molecules would be deposited in $a_4$ configurations without interaction between them. This configuration appears to be more stable than those corresponding to the deposition of a T-shaped or SE dimer. The negative sign obtained for the $a_4-cd$ structure shows that the latter is the most stable investigated configuration as it maximizes the interaction between each benzene molecule and the graphene surface while maintaining some stabilizing intermolecular interactions between the two benzene units. Again, this conclusion holds at both DFT and DFTB levels. It should also be noted that the values of the interaction energies are of similar order for these two levels of calculation.

6. Conclusions

In the present paper, we have reported a new implementation of periodic boundary conditions in the DFTB code deMonNano, as only the $\mathrm{\Gamma }$-point approximation was available in the previous version of the code. An originality of our scheme is the inclusion of atomic charge corrections, which improves the description of intermolecular coulomb interactions. It allows the recovery of a reasonable description of molecular clusters, as shown in the particular case of benzene dimers in this work. Dispersion corrections are also mandatory for a proper description of such interactions, and we have benchmarked two empirical correction schemes. One of them produces the best C-C bond distance in graphene and graphite, whereas the second one provides the best interlayer distance in graphite according to previous reference calculations, experiments and new DFT calculations performed with the DFT-D3 dispersion-corrected functional.

Benzene monomer and dimers have been optimized at the DFTB and DFT levels, providing the following similar trends. For the deposition of a single benzene monomer on a graphene sheet, the adsorption of benzene centered on top of graphene carbon atom or C-C bond leads to almost degenerated structures, which are by far more stable than the superimposition of the benzene on top of a graphene aromatic cycle. The most stable one at the DFTB level has been selected to build initial conditions for benzene dimers deposition on graphene. The structural energy gap between the most stable T-shaped dimer and less stable Sandwich-like dimer is divided by two when the cluster is supported on graphene. The supported Parallel-Displaced structure appeared to be unstable at the DFTB level, resulting in a structure where the two benzene are deposited close to each other on the graphene surface. This structure is the most stable one of our calculations at DFT and DFTB levels and also more stable than the deposition of two benzene monomer at infinite distance, which is not the case of the deposited sandwich or T-shaped dimers.

As a conclusion, we have shown the ability of the new implementation to characterize properties of molecular clusters deposited on surfaces, opening the path to dynamical simulations that allow probing the role of deposition conditions (e.g., cluster/support collision energy) and temperature on the structure and stability. This work validates the implementation of the improved sampling of the reciprocal space (k-points) in deMonNano, which will allow to perform very precise calculations, as needed to study properties for which a very fine structure of the Brillouin zone is required, such as electronic transport, magnetism or topological states of matter [57].