Modeling Interactions between Graphene and Heterogeneous Molecules

: The Lennard–Jones potential and a continuum approach can be used to successfully model interactions between various regular shaped molecules and nanostructures. For single atomic species molecules, the interaction can be approximated by assuming a uniform distribution of atoms over surfaces or volumes, which gives rise to a constant atomic density either over or throughout the molecule. However, for heterogeneous molecules, which comprise more than one type of atoms, the situation is more complicated. Thus far, two extended modeling approaches have been considered for heterogeneous molecules, namely a multi-surface semi-continuous model and a fully continuous model with average smearing of atomic contribution. In this paper, we propose yet another modeling approach using a single continuous surface, but replacing the atomic density and attractive and repulsive constants in the Lennard–Jones potential with functions, which depend on the heterogeneity across the molecules, and the new model is applied to study the adsorption of coronene onto a graphene sheet. Comparison of results is made between the new model and two other existing approaches as well as molecular dynamics simulations performed using the LAMMPS molecular dynamics simulator. We ﬁnd that the new approach is superior to the other continuum models and provides excellent agreement with molecular dynamics simulations.


Introduction
The physical adsorption of various molecules onto substrates, a process known as physisorption, has been shown to be a viable mechanism in gas storage [1,2], pollutants capture [3,4] and binding of biological molecules including proteins [5] and DNA [6]. Such adsorption works primarily through the van der Waals interactions and does not involve chemical binding, as opposed to chemisorption. Carbon allotropes are often considered as candidates for adsorbents in physisorption for molecular capture or storage [7,8]. Apart from experimental studies, techniques used to model physisorption include molecular dynamics simulations [9][10][11][12], Monte Carlo methods [11,13,14] and continuum mathematical models [15][16][17][18] to name only a few. In particular, mathematical modeling of the interactions between various molecules and carbon surfaces has been extensive [19][20][21][22][23][24][25][26]. The standard approach for this modeling is to use the Lennard-Jones potential together with a continuum approximation [27], which assumes a uniform distribution of atoms over the surface of the molecule or throughout the volume of the molecule. The use of the Lennard-Jones potential is justified since the dominant forces present within the interactions are van der Waals forces. This approach allows the potential energy between the two interacting molecules to be expressed as a double integral over two surfaces or throughout two volumes of the molecules [28]. This integral can usually be evaluated in terms of special functions, giving rise to a closed-form analytical expression, which is useful for the rapid determination of numerical results for predictive purposes, such as determining the size of nanotubes that are capable of encapsulating other molecular structures [27]. The surface approximation of the molecules plays a key role in determining the accuracy of the resulting expression [27], which makes this modeling technique difficult when considering molecules that are irregularly shaped or have a non-uniform distribution of atoms.
The current approach to overcoming this difficulty in modeling complex molecules, such as DNA or proteins, interacting with other molecules is to approximate their surfaces or volumes and assume that the atoms are homogeneously distributed throughout the structure [22,24,29,30]. This allows for the continuing use of constants for atomic density and the attractive and repulsive coefficients. A potential problem with this approach is the homogeneous smearing of atoms which may produce inaccurate results when the molecule being approximated is heterogeneous (i.e., comprising different atomic species).
In this paper, we use interaction functions recently proposed by Stevens et al. [31] to replace the Lennard-Jones constants (A and B). In addition, this paper develops a density function to replace the constant atomic surface (or volume) density in the continuum approximation. This technique greatly improves the accuracy for the energy calculation since it takes into account the heterogeneity of the molecules. To demonstrate the use of this technique, we consider the interaction of a coronene molecule with a graphene sheet. Due to the symmetry of a coronene structure (see Figure 1), we assume the attractive and repulsive functions (A(r) and B(r)) to depend on the radial distance r of the coronene molecules. We then compare the results of the interaction energy with three different continuum models. Model (i) involves a homogeneous smearing of atoms with constant atomic density, attractive and repulsive constants. Model (ii) assumes a constant atomic density but with attractive and repulsive functions which are functions of r. Finally, Model (iii) is a model with atomic density, attractive and repulsive as functions of r. These results are also benchmarked with molecular dynamics simulations performed using LAMPPS. In the following section, we describe the basic model including the Lennard-Jones potential, the formulation of the integral including the density and interaction functions and how it applies to the interaction of a coronene with a sheet of graphene. In the section thereafter, we detail the results for the various approaches by comparing the results for the three models with data obtained from molecular dynamics simulations. Some brief conclusions are presented in the final section of the paper. Appendices A and B include some technical details that are required in the determination of formulae for the disk and plane integral described in the methods.

Methods
Using the Lennard-Jones potential, the interaction energy between two molecules can be expressed as a double surface integral over the surfaces the interacting molecules, namely where η 1 and η 2 are the atomic surface densities of the two molecules, A and B are the attractive and repulsive constants, respectively, and ρ is the typical distance between the surface elements dS 1 and dS 2 . Equation (1) is referred to as a continuum approach, which is based on pair-wise summation of the interatomic interactions between both molecules [28]. In general, this approach assumes that atoms are uniformly distributed over the surface of the molecules, and thus the atomic surface density can be obtained from η = N/S, where N is the total number of atoms on the molecule and S is its surface area. Similarly, the attractive and repulsive constants A and B are calculated based on the types of interacting atoms on the molecules [27].
The continuum approach has been used to investigate the interactions between many types of molecules, especially those with regular shaped structures. While this approach works well with homogeneous materials which comprise same types of atoms, it requires significant improvement and modification to give a more accurate result for heterogeneous structures comprising more than one types of atoms.
Here, we propose an improved continuous model for a heterogeneous molecule by replacing constants A and B in (1) with interaction functions, Here, r is the vector of variables that parameterize the heterogeneity of the interaction between the molecules. In the case of a heterogeneous molecule interacting with a homogeneous molecule, the heterogeneity is parameterized only in terms of the heterogeneous molecule. Furthermore, it is assumed that both A (r) and B (r) have the same functional form. This allows for the use of the standard computational technique, solving some K n for all integer n where Thus, the expression in (2) can be written using (3) as E = η 1 η 2 (−K 3 + K 6 ). For the coronene-graphene interaction, we consider graphene as a homogeneous plane of carbon atoms. For coronene, we model as a disk for which carbon and hydrogen atoms dominate the middle region and its perimeter, respectively. As the coronene is radially symmetric, we assume the interaction functions vary only along the radial direction, r.
Due to the structure of coronene as shown in Figure 1b and that the interaction of coronene with graphene is made up of carbon-carbon and carbon-hydrogen interactions, we assume that the interaction functions satisfy the conditions: where the Lennard-Jones constants for carbon-carbon and carbon-hydrogen interactions are as given in Table 1 and R is the radius of the coronene molecule, which is taken to be 4.79 Å [32]. In this paper, we adopt the interaction functions to be in the forms of sigmoidal arctan functions as we assume that the interaction varies similar to a step-function but with smooth transition across the radial direction of the molecule. Specifically, we have where α and β are determined using Equation (4). We choose a large value for m in order to mimic the step-function profile, and the value of r 0 is chosen to be located between the connection of final carbon atom and the hydrogen atom at the perimeter. Thus, we use m = 500 and r 0 ∈ (0.785, 1). The profile of the interaction functions is shown in Figure 2a.

Atomic Density Function
The atomic densities of both surfaces are typically assumed to be uniform, meaning that η 1 and η 2 are constants and are not included in the integrals, as shown in (1). Here, we assume that the atomic surface density of the coronene molecule is not a constant, i.e., replacing η 1 with a function η 1 (r), thus bringing it into the integral along with the interaction functions, Due to the radial symmetry of the coronene molecule, we also express the atomic density as a function of the radial distance, r. Using the atomic densities calculated in the semi-continuous ring model developed by Tran-Duc et al. [32], a quadratic interpolating polynomial is fitted and used as the atomic density function in our calculation. These densities are presented in Table 2 and the polynomial is plotted in Figure 2b.
To solve the energy integral for a general case, we use to express the atomic density function, where N here is the degree of the interpolating polynomial.

Disk and Plane Interaction
Coronene comprises 24 carbon atoms and 12 hydrogen atoms, arranged in a circular honeycomb structure with carbon atoms in the central region and the hydrogen atoms at the perimeter. Coronene has previously been modeled as concentric rings of carbon and hydrogen atoms [32]. Here, we represent a coronene molecule by a circular disk of radius a centered on the z-axis, and tilted by some angle φ ∈ [0, π /2]. We also assume that the centre of the disk is located at a distance δ > a sin φ away from a sheet of graphene which is assumed to lie flat on the xy−plane.
Next, we define the integral K n as where ρ 2 = (x − ar cos φ cos θ) 2 + (y − ar sin θ) 2 + (δ − ar sin φ cos θ) 2 . Evaluating (8), we find that the integral is largely the same as a circular ring interacting with a plane, as the integral in r is independent of the other variables. Thus, we refer the reader to the work of Tran-Duc et al. [32] for detailed computation of most of the K n integral. The expression for Equation (8) prior to solving the integral in r is where Detailed derivation of Equation (10) can be found in Appendix A.

Results and Discussion
Here, we compare the energies obtained from the three different models: (i) the fully homogeneous disk with constant interaction and density coefficients; (ii) the heterogeneous disk with functional interaction and constant density coefficients; and (iii) the heterogeneous disk with functional interaction and density coefficients, against the results from molecular dynamics (MD) simulations. The MD simulations are performed using the Lennard-Jones potential with a cut-off distance of 14 Å and the graphene sheet used has radius over 20 Å larger than the radius of the coronene molecule (4.8 Å), which justifies the assumption of an infinite plane for the sheet.
Fixing the tilt angle at φ = 0 implies that the coronene molecule lies parallel to the graphene sheet. As seen from Figure 3, the differences in the energy profiles obtained from the three models are obvious. Here, Model (i) agrees excellently with the MD simulation, whereas Models (ii) and (iii) using interaction functions have discrepancies. We find that this is due to Model (i) and the atomistic case being equal when the coronene is parallel to the graphene sheet, which is shown in Appendix B.
Energy, E (kcal/mol) Where Model (i) begins to fail appears to depend on the distance between the coronene and graphene. Figure 4a shows that it fails at φ ≥ 3π /64 when δ =3.5 Å and Figure 4b shows that it fails at φ ≥ 3π /32 when δ =5 Å. Conversely, we notice that Models (ii) and (iii) appear to maintain or increase their accuracy at these fixed distances as φ increases. For φ > 0, we see that Model (i) becomes very inaccurate for δ < δ min + 0.5 where δ min is the distance at which the energy is minimized (E min ). This can be seen in Figure 5 for various values of φ. Model (ii) produces energy curves that are clearly an improvement over the homogeneous model, as shown in Figure 5a,b. The minimum energy distances presented in Table 3 show this to be the case as the model is consistently within 0.1 Å of the MD simulation. This is in contrast to the error between the simulation prediction and Model (i), which is initially small at 0.01 Å when φ = 0 but very quickly increases to 0.29 Å by φ = π /2. Figure 5. Energy profiles for coronene at fixed tilt angles: (a) φ = π /6; (b) φ = π /4; (c) φ = π /3; and (d) φ = π /2. The combination of the interaction functions and density function seems to correct the discrepancy in energy estimation. This is shown prominently in Figure 5c,d as Model (iii) follows Model (ii) closely until the distance nears δ = δ min . On calculating the E min error, E = |E min,model − E min,simulation | , at each tilt angle φ, we find that this error for the homogeneous Model (i) decreases from E = 2.26 kcal mol −1 at φ = π /6 to E = 1.37 kcal mol −1 at φ = π /2. An improved error is obtained for Model (ii), for which E = 1.52 kcal mol −1 at φ = π /6 to E = 0.78 kcal mol −1 at φ = π /2. However, Model (iii) gives the smallest error over this range of tilt angles with E = 0.54 kcal mol −1 at φ = π /6 to E = 0.17 kcal mol −1 at φ = π /2. The error values in δ min and E min for the three models indicate that including an interaction function (Model (ii)) increases the accuracy over the purely homogeneous model (Model (i)), which is improved upon further by incorporating a density function (Model (iii)). Thus, unless the coronene molecule alignment is perfectly parallel to the graphene sheet, the homogeneous model of coronene should not be adopted.

Conclusions
While the models presented here do not perfectly match the MD results for the energy of the coronene-graphene interaction, it is clear that the fully homogeneous Model (i) gives the largest error in most configurations. Our proposed method could lead to a study to improve the homogeneous models used for the interactions involving a much more complex molecular structures than the coronene molecule. The introduction of simple functions to approximate the interaction strengths and the atomic densities will likely lead to improved energy estimations, which will be the subject of future research.

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

Abbreviations
The following abbreviations are used in this manuscript:

LJ
Lennard-Jones MD Molecular dynamics Appendix A. Evaluation of the Integral I Defined by (10) To evaluate the integral (10), we use the interaction function f (r) = α arctan (m (r 0 − r)) + β. We need to compute The second part is easily evaluated since I 2 = 1 0 r k dr = 1 /(k+1). For the first part, we make the substitution u = m (r 0 − r) in order to make the argument of the arctan function easier to manage, and now we use the binomial theorem to split up the polynomial term in u into the series, From [33], we have, and the second integral has two possible forms depending on whether m is even or odd. For m odd, we obtain while, for m even, we obtain To make use of these forms, we split the finite sum in Equation (A2) into two sums of odd and even powers, so that, on applying the integral (A3) and noting that there are no singularities over the interval, we obtain We then apply (A4) to deduce Although this looks complicated, due to the series being finite, the expression is easily computable by any mathematical software package, such as MAPLE™ or Mathematica™. and since this is a finite sum of the above integrals based on the degree of the polynomial approximating the density, the solution to this is the same as in Equation (A5).

Appendix B. Explanation to Why MODEL (i) Is Accurate at φ = 0
In this appendix we show that the energy expression for Model (i) matches the energy expression given by the atomistic model of coronene solely when the tilt angle is zero.
Appendix B.1. Model (i) N = N C + N H is the total number of atoms in coronene with N C the number of carbon atoms and N H the number of hydrogen. η coro = N SA coro is the atomic density of coronene and SA coro is the surface area of the approximated coronene surface. Thus, the interaction energy between a coronene molecule and a parallel graphene sheet is computed by

Appendix B.2. Atomistic Model
The interaction energy between a coronene molecule and a parallel graphene sheet is the sum of the interaction energies between each atom and the graphene sheet, where ρ k is the distance between the kth carbon or hydrogen atom and the graphene sheet.
As with the homogeneous case, at φ = 0, ρ k = ρ ∀k. Thus, E coro-graph = η graph which is equal to the expression derived from the homogeneous approximation.