AIRBED: A Simpliﬁed Density Functional Theory Model for Physisorption on Surfaces

Dispersion interactions are commonly included in density functional theory (DFT) calculations through the addition of an empirical correction. In this study, a modiﬁcation is made to the damping function in DFT-D2 calculations, to describe repulsion at small internuclear distances. The resulting Atomic Interactions Represented By Empirical Dispersion (AIRBED) approach is used to model the physisorption of molecules on surfaces such as graphene and hexagonal boron nitride, where the constituent atoms of the surface are no longer required to be included explicitly in the density functional theory calculation but are represented by a point charge to capture electrostatic eﬀects. It is shown that this model can reproduce the structures predicted by full DFT-D2 calculations to a high degree of accuracy. The signiﬁcant reduction in computational cost allows much larger systems to be studied, including molecular arrays on surfaces and sandwich complexes involving organic molecules between two surface layers.


Introduction
The study of 2-dimensional (2D) arrays of organic molecules on surfaces is an increasingly important area of research owing to their potential applications in a range of fields including electronic and optoelectronic devices. [1][2][3][4] Understanding the interactions that underpin the organisation of these arrays is key to realising their potential. These interactions can include hydrogen bonding, covalent bonding, dispersion forces and metal coordination, which may all have a significant role in the structure of the array. [5][6][7][8] One area of particular focus is the study of organic arrays on surfaces such as graphene and hexagonal boron nitride (hBN). These arrays can be imaged using scanning tunnelling microscopy and atomic force microscopy, which can provide high-resolution images of 2D molecular and supramolecular organization under ultra-high-vacuum conditions and at atmospheric pressure. [9][10][11][12][13] When adsorbed on insulating surfaces, the fluorescence of the adsorbed molecular layers can also be measured providing a direct link to the optoelectronic properties of the arrays. [14][15][16][17][18][19][20] Recent experimental work has shown shifts in the fluorescence emission bands arising from adsorption on hBN can be measured, and it has also been demonstrated that these shifts can arise from changes in the molecular structure that occur on absorption, and other factors such as screening effects that depend on the refractive index of the substrate. [14][15][16][17][18][19][20] Quantum chemical calculations, for example density functional theory (DFT), can make a significant contribution to characterising the structure of adsorbed molecular arrays and also their optoelectronic properties. An inherent limitation of molecular quantum chemical calculations for these systems is the computational cost of treating a large organic molecule and the model of the surface so that only a single molecule may be typically considered. 18 Clearly it is desirable that many adsorbed organic molecules are included in the calculation to capture the effects of the larger molecular array, and this motivates the search for computationally less expensive approaches that include the important interactions that can be used to model the structure of these systems.
Quantum chemical methods, including coupled cluster theory, have been used to determine the absorption energies of small molecules on molecular models of the graphene and fluorographene surfaces. [21][22][23] A symmetry adapted perturbation theory analysis showed that the dominant contribution to the adsorption energy arises from dispersion 21 which is consistent with other work which also found that dispersion is the most significant interaction between the insulating surface and adsorbed array for systems where the molecular arrays are physisorbed on the surface. 18 The next largest contribution to the adsorption energy came from the electrostatic interactions. In this work we focus on DFT since this can be most readily applied to larger systems. Incorporating dispersion forces within DFT is an active area of research. [24][25][26][27][28][29][30][31][32][33] However, the most computationally efficient approaches are empirical dispersion models. An example of this approach is the DFT-D2 dispersion correction of Grimme,28 in this approach the total energy is expressed as where E KS−DF T is the usual DFT energy according to the chosen functional, and E DISP is the dispersion energy that is typically given by Here R AB and R 0 AB are the internuclear separation and sum of the van der Waals radii of atoms A and B respectively, C AB 6 is the dispersion coefficient for atom pair AB, s 6 is a scaling factor and there are N atoms in the system. The role of the damping function, f dmp (R AB ) is to avoid double counting of electron correlation effects and the near-singularites as R AB → 0.
In this paper, we describe a modification to the damping function so that it also describes the Pauli repulsion between the atom pairs, with the contribution from electrostatic interactions included by assigning a partial charge to the atomic sites of the surface. The resulting model is then used to study the structure of molecular arrays on the graphene and hBN surfaces, and organic molecules trapped between two layers of graphene of hBN. 34

Computational Methods
In the approach used here, the atoms of the molecular environment, in this case the surface atoms, are not included in the DFT calculation but are represented by points in space that can be assigned a charge q. The interaction between the surface and the molecule is described by a modification of the dispersion correction to include repulsion. This is achieved through a modification of the damping function; specifically the surface-molecule interaction E DISP is replaced by E vdW , where for the N s surface "atoms" and N m molecule atoms. In this model the repulsion at short range is described by an exponential function. This functional form for the repulsion is consistent with the work of Wheatley and Price that that recognised the proportionality between the exchange-repulsion energy and the charge density overlap. 35 The R 0 AB values used for this contribution to the energy are derived from the experimental van der Waals radii, 36 which tend to be greater than the values used in the standard DFT-D2 correction.
The difference in atomic radii used for the two methods reflects the fact that in the conventional approach, the DFT calculation includes interactions with substrate atoms, which will contribute to the intermolecular interaction, and this interaction is no longer present.
This has also been observed in other approaches that combine empirical potentials with the Grimme-type dispersion correction. 37 A value of d=20.0 is used, which is unchanged from the original damping function and s 6 and the C AB 6 coefficients are also unchanged. An additional parameter α is introduced which allows some additional flexibility to tune the surface-molecule interaction. This parameter has values of 0 and 1.2 for graphene and hBN, respectively.
We emphasise that for the interaction between the atoms of the adsorbed molecule(s), the original, unmodified dispersion correction is used. A comparison between the original dispersion contribution, E DISP , and the modified function, E vdW , is shown in Figure 1 for approach where the interaction with the surface is treated by an empirical potential has been used previously, 38 however the approach presented here can be viewed as a QM/MM approach that is integrated within the DFT calculation. This feature provides the approach with the potential to be applied to a wide range of systems that would normally be studied with DFT with dispersion.
The AIRBED model is assessed through a comparison with the structure of a single molecule adsorbed on the graphene and hBN surfaces. A set of four molecules; 4,4- The DFT-D2 calculations predict that the molecules are further from the surface of graphene compared with hBN, and this trend is replicated by the AIRBED calculations.
The average error in the predicted heights from AIRBED is less than 0.1Å compared with the DFT-D2 calculations, with the largest discrepancy of 0.18Å occurring for PTCDI on graphene. As expected, the RMSD for the rigid planar molecules PTCDI and H 2 Pc are very small and the resulting structures are indistinguishable. For the more, flexible molecules 4,4-DIB and TIPB there is an increase in the RMSD. However, an RMSD of 0.48Å still represents a close match between the structures. This is illustrated in Figure 3 which shows the DFT-D2 and AIRBED optimised structures overlaid for the two systems with largest RMSD. There are no striking visual differences but a closer inspection shows some misalignment of the phenyl rings.   An extension of these systems which may be invisaged are so-called "sandwich" complexes where the organic layer is positioned between two layers. Within a full DFT-D2 treatment, the number of atoms required to model both layers would make the calculation very computationally demanding. However, within the AIRBED approach the inclusion of the layer does not add any practical increase in the cost of the calculation. Figure 7 shows the energy of the PTCDI molecule within two layers of graphene as the interlayer separation is varied. ference between the S 0 and S 1 states for the two surface models agree to within 0.01 eV.
This indicates that the findings for the validity of the model for ground state structures also applies for excited states.

Conclusions
A modification of the damping function used in the standard DFT-D2 method to capture repulsion in addition to dispersion interaction has been introduced. This provides a QM/MM approach integrated within the DFT calculation. The resulting AIRBED model is shown to reproduce closely the structures of molecules physisorbed on graphene and hBN surfaces predicted by full DFT-D2 calculations at a vastly reduced computational cost. The model has been used to study molecular arrays adsorbed on hBN and sandwich complexes. These systems are of interest owing to their potential uses in optoelectronic devices, however, the AIRBED approach can be applied to a much wider array of systems including the study of the molecular dynamics of the adsorbed molecules, and these are the subject on on-going investigations.

Conflicts of interest
There are no conflicts of interest to declare.