Interaction between Graphene-Based Materials and Small Ag, Cu, and CuO Clusters: A Molecular Dynamics Study

The excessive use of antibiotics has contributed to the rise in antibiotic-resistant bacteria, and thus, new antibacterial compounds must be developed. Composite materials based on graphene and its derivatives doped with metallic and metallic oxide nanoparticles, particularly Ag, Cu, and Cu oxides, hold great promise. These materials are often modified with polyethylene glycol (PEG) to improve their pharmacokinetic behavior and their solubility in biological media. In this work, we performed molecular dynamics (MD) simulations to study the interaction between small Ag, Cu, and CuO clusters and several graphene-based materials. These materials include pristine graphene (PG) and pristine graphene nanoplatelets (PGN) as well as PEGylated graphene oxide (GO_PEG) and PEGylated graphene oxide nanoplatelets (GO-PEG_N). We calculated the adsorption energies, mean equilibrium distances between the nanoparticles and graphene surfaces, and mean square displacement (MSD) of the nanoclusters. The results show that PEGylation favors the adsorption of the clusters on the graphene surfaces, causing an increase in adsorption energies and a decrease in both distances and MSD values. The strengthening of the interaction could be crucial to obtain effective antibacterial compounds.

Molecular dynamics (MD) simulations allow the study of materials from a molecular viewpoint, providing information on an atomic scale usually inaccessible to experimental techniques [33]. The interaction between graphene-based materials and different nanoparticles and polymers has been described in many research works by MD [34][35][36][37][38][39][40][41][42]. Several

Models and Calculation Method
MD simulations were done with the Forcite module of the Materials Studio 9 software [49] in the NVT ensemble (constant number of particles N, constant volume V, and constant temperature T) at 298 K during 10 ns of simulation length. The interactions between atoms were calculated using the COMPASS forcefield, which is a force field parameterized using experimental and ab initio results that enables an accurate prediction of various gas-phase and condensed-phase properties of several organic and inorganic materials [50], including graphene and graphene-based materials [44,51,52]. The temperature was kept constant using a Nosé-Hoover thermostat . [53,54]. Electrostatic and van der Waals interactions were calculated using an atom-based summation method with a cut-off of 12.5 Å.
The starting configurations used for the production trajectories are shown in Figure 1 (graphene monolayers) and Figure 2 (nanoplatelets). A basal graphene sheet with a width of 46.5 Å and a length of 56.2 Å was used for all structures. To build these starting structures, all the systems were subjected to an annealing procedure composed of NVT MD simulations and further steepest descent and conjugate gradient energy minimization cycles at each temperature. This annealing procedure allows exploring the conformational space for low-energy structures, which will be used as starting configurations, by periodically increasing and then decreasing the temperature (from 300 up to 1000 K) of a classical dynamics trajectory, to avoid trapping the structure in local energy minima. For the sake of comparison, PG, PGN, GO_PEG, and GO-PEG_N were constructed. The PEGylated structures were functionalized with eight short chains of PEG (degree of polymerization n = 10). In the interest of saving computational time, only two layers of graphene were used to model the nanoplatelets. In the PEGylated nanoplatelets, both layers had the same grafting density, i.e., the number of PEG chains attached to the edges of the individual GO layers. In these starting structures, the mean final distances between the basal planes of graphene were 3.7 and 6.9 Å for PGN and GO-PEG_N, respectively. The clusters used to represent the metal nanoparticles were composed of thirteen atoms arranged in an icosahedral shape. This number was chosen because of the well-known stability of this regular icosahedron geometry for Cu and Ag nanoclusters [55,56]. For CuO nanoparticles, a Cu 6 O 6 cluster consisting of two hexagonal layers, one above the other, was chosen. This structure was found to be stable by a density functional theory calculation [57]. No constraints were imposed on the systems, so that all molecules could move freely over the whole simulation length. At the beginning of the calculation, the clusters were placed far from the surfaces, at distances of about 8 Å, as can be seen in the lower part of both Figures, and they were left to evolve with time.
clusters were placed far from the surfaces, at distances of about 8 Å, as can be seen in the lower part of both Figures, and they were left to evolve with time.       The adsorption energy of the clusters on the surface was calculated from the following equation: where E GBM+cluster is the mean equilibrium potential energy of the graphene-based material (GBM = graphene-based materials = monolayers or nanoplatelets) interacting with the cluster, and E GBM and E cluster are the mean equilibrium potential energies of the isolated graphene-based materials and clusters, respectively (all the energies were calculated as average values from three independent 10 ns simulations. In turn, the result of each individual simulation is the average of the energy values obtained during the last 5 ns of simulation). As the energy values are negative, the more negative the adsorption energy is, the stronger the interaction. Previous studies have shown that this adsorption energy could be a good estimate of the binding strengths of adsorbate-adsorbent systems [58]. The mean square displacement (MSD) of the clusters, which measures the spatial extent of clusters random motion, was calculated from their positions over simulation length according to: where r i (0) is the reference position of the cluster, r i (t) is its position at time t, and N is the number of atoms of the cluster. Thus, MSD was calculated from the particle positions obtained from the trajectory using a time step of 5 ps. MSD can be related to time through the following equation: where K α is the generalized diffusion coefficient and α is the diffusion exponent [59,60]. The value of this exponent is related to the mechanism of diffusion.

Adsorption Energies
The adsorption energy values calculated from Equation (1) are presented in Table 1. From the values listed in Table 1, it can be concluded that the strength of the interaction depends on both the chemical nature of the cluster and the functionalization of the graphene-based material. All the clusters interact more weakly with pristine materials than with functionalized materials in both monolayers and nanoplatelets. The Cu 6 O 6 cluster is the most attracted by the surface and the Cu 13 cluster is the least attracted. Ag 13 lies in between. The introduction of PEG chains favors the adsorption process, increasing the adsorption energy. In addition, the nanoplatelets seem to increase the adsorption energy Here, we must note the qualitative nature of these results. It is well known that cluster geometry and size influence their adsorption energy and reactivity on graphene surfaces, and several authors have shown that quantum mechanical calculations are needed to fully understand and quantify this kind of interaction [61]. Even when using density functional theory calculations, the results vary widely (adsorption energies ranging from a few kcal/mol to a hundred kcal/mol and cluster-graphene surface distances ranging from 2 to 3 Å) depending on the molecular models, the calculation method, or the number of atoms included in the nanocluster [62][63][64]. What is clear is that metal d-orbitals are involved in charge transfer processes between the clusters and the graphene surface, although this charge transfer is weak in some cases, and the bonding can be mainly attributed to van der Waals forces. Cu and Ag are examples of that, as it has been demonstrated by numerous calculations [65,66]. On the other hand, experimental results are diverse and show that the structure of the material after the interaction depends on many factors, including the initial graphene form and the methods used for its synthesis [67]. In addition, due to limited computational resources, our simulations were performed in vacuum, and we are aware that this may influence the results. Thus, the small clusters and the graphene model used in this work represent a first step to compare qualitatively the three types of particles, and we realize that they cannot capture the whole complexity of real systems. Cu 6 O 6 interacts strongly with both pristine and functionalized surfaces, although the adsorption process is improved by the presence of PEG chains, as already observed for Cu and Ag clusters. We think the atomic charges of the cluster increase the electrostatic attraction between the cluster and surface, although the exact nature of the interaction should be verified by more accurate quantum mechanical calculations. Sun et al. [68] determined that during the process of nucleation of graphene on Cu-based substrates, the adsorption energies of C were larger for copper oxides than for metallic Cu. Ko et al. [69] prepared CuO nanoparticles covered with a monolayer graphene shell, in which strong C-O-Cu links were formed. Using density functional theory calculations, Mohammadi-Manesh et al. [70] calculated the binding energy of different configurations of adsorption of Cu and CuO on graphene surfaces and found that the interaction between Cu and graphene is physical and the interaction between CuO and graphene is chemical.
Even if the calculation method used in this work cannot describe the formation of chemical bonds between the cluster and the graphene surfaces, the results obtained agree with those obtained from more accurate quantum chemical calculations.

Cluster-Surface Distances
The mean equilibrium distances between the cluster and nearest surface atoms are shown in Table 2. The reduction of the distances used in the starting structures (8 Å) indicates that the clusters approached the graphene surface as the simulation progressed and attained an equilibrium value toward the end of the calculation. The values listed in Table 2 show that the mean distance between the clusters and the pristine surfaces is about 2.9 Å for both monolayers and nanoplatelets. The introduction of PEG chains brings the cluster closer to the graphene surface leading to reduced distances of about 2.4-2.5 Å. Cluster-pristine material distances are larger than cluster-functionalized material distances, which correlates well with improved adsorption energies for GO_PEG-containing materials.
The distances shown in Table 2 for the pristine systems are similar to those found in previous quantum mechanical [65,66] and MD [71,72] simulations of the interaction between graphene and Cu and Ag. DFT results give distances of 2.3-3.9 and 2.8 Å for Cu and Ag, respectively. The distance found by MD between Ag and graphene is 3 Å. MD results for Cu show a range of binding distances varying between 2 and 3.4 Å, depending on the orientation of Cu relative to graphene. They are related to weak bonding (physical adsorption) attributed to the full occupancy of metal d orbitals [73,74]. The radial distribution functions (RDF) of clusters and graphene-based surfaces are shown in Figure 3. These plots give us an idea of the spread of distances between the clusters and surfaces. All plots present maxima between 2.5 and 5-6 Å, and there are no large differences between the different models. It must be noted that the mean equilibrium distances shown in Table 2 represent distances between the cluster and nearest surface atoms, which must not necessarily coincide with maxima in Figure 3. In the case of the Cu 6 O 6 cluster, the plots corresponding to the Cu atom (shown in red) are slight shifted to larger values, which is in accordance with the results shown in Table 2. Figure 3 also shows that in the PEGylated systems, distances are somewhat smaller.  Figure 4) and the nanoplatelet surface ( Figure 5). Only the result obtained from one of the simulations is depicted here. In the case of the GO_PEG surfaces, it was observed that the clusters could stick to any lateral PEG chain of the layers, and once stuck, they remained bonded to them at a distance of about 2.5-2.6 Å in all cases (see Table 2). In   Figure 4) and the nanoplatelet surface ( Figure 5). Only the result obtained from one of the simulations is depicted here. In the case of the GO_PEG surfaces, it was observed that the clusters could stick to any lateral PEG chain of the layers, and once stuck, they remained bonded to them at a distance of about 2.5-2.6 Å in all cases (see Table 2). In the case of CuO clusters, this distance corresponds to O-PEG distances. In all cases, the Cu atom in the cluster was located at somewhat larger distances. PEG is frequently used to prepare metal and metal oxide nanoparticles in solution [75][76][77], as it is well known that PEG molecules strongly adsorb on metal nanoparticles surfaces by coordination through the ether bond of PEG and prevent their aggregation [78]. The capping capability of PEG could explain the strong attraction it exerts on the clusters and their subsequent immobilization on the graphene surfaces.

Mean Square Displacement
Although a detailed analysis of the mechanisms of cluster diffusion is beyond the scope of this work, we think that the MSD of the cluster center of mass can help shed

Mean Square Displacement
Although a detailed analysis of the mechanisms of cluster diffusion is beyond the scope of this work, we think that the MSD of the cluster center of mass can help shed more light on the differences observed in adsorption energies. Total MSD is plotted versus simulation length in Figures 6 and 7 for monolayers and nanoplatelets, respec-

Mean Square Displacement
Although a detailed analysis of the mechanisms of cluster diffusion is beyond the scope of this work, we think that the MSD of the cluster center of mass can help shed more light on the differences observed in adsorption energies. Total MSD is plotted versus simulation length in Figures 6 and 7 for monolayers and nanoplatelets, respectively. Larger MSD values indicate greater cluster mobility.  Cu13 and Ag13 clusters interacting with the pristine materials present the highest Cu 13 and Ag 13 clusters interacting with the pristine materials present the highest MSD values. Both PG and PGN weakly attract both clusters, which keep constantly moving over the material surface. However, the functionalization of the graphene layers drastically reduces their mobility (see bottom plots in Figures 6 and 7. The PEG chains act by trapping the clusters and blocking their displacement over the surface. The Cu 6 O 6 cluster shows a completely different behavior as its MSD values are low (blue line plots) in both the pristine and functionalized materials.
To show more clearly the cluster movement, the diffusion trajectories of the cluster center of mass over the graphene surface during the last 250 ps of simulation length for both PG and PGN are displayed in Figures 8 and 9, respectively. Some video representative of the different cluster dynamics is included in the Supplementary Materials (Videos S1, S2, and S3).
The MSD values in all directions are shown in Figures 10 and 11 for monolayers and nanoplatelets, respectively. The graphene surfaces are oriented parallel to the xy plane. As can be inferred from these figures, the MSD in the z direction is approximately zero (green line), as once the cluster approaches the basal plane, no vertical displacement from the surface is observed. MSD xx and MSD yy components have non-zero values due to the movement of the cluster over the surface and parallel to it during the whole simulation length. In the functionalized structures, the cluster may carry out slight movements between or along PEG chains (see Figure 12), although it remains bonded to them. Therefore, MSD xx and MSD yy present low non-zero values. An example of this dynamic behavior is included in the Supplementary Material (Video S4).   moving over the material surface. However, the functionalization of the graphene layers drastically reduces their mobility (see bottom plots in Figures 6 and 7). The PEG chains act by trapping the clusters and blocking their displacement over the surface. The Cu6O6 cluster shows a completely different behavior as its MSD values are low (blue line plots) in both the pristine and functionalized materials.
To show more clearly the cluster movement, the diffusion trajectories of the cluster center of mass over the graphene surface during the last 250 ps of simulation length for both PG and PGN are displayed in Figures 8 and 9, respectively. Some video representative of the different cluster dynamics is included in the Supplementary Material (Videos S1, S2, and S3).     The MSD values in all directions are shown in Figures 10 and 11 for monolayers and nanoplatelets, respectively. The graphene surfaces are oriented parallel to the xy plane. As can be inferred from these figures, the MSD in the z direction is approximately zero (green line), as once the cluster approaches the basal plane, no vertical displacement from the surface is observed. MSDxx and MSDyy components have non-zero values due to the movement of the cluster over the surface and parallel to it during the whole simulation length. In the functionalized structures, the cluster may carry out slight movements between or along PEG chains (see Figure 12), although it remains bonded to them. Therefore, MSDxx and MSDyy present low non-zero values. An example of this dynamic behavior is included in the Supplementary Material (Video S4).   The MSD values corroborate the trends followed by the adsorption energies: the more dynamic the cluster is, the lower the adsorption energy. In agreement with our results, Gervilla et al. [79,80] studied the diffusion of Cu and Ag adatoms and clusters on graphene, using ab initio and classical molecular dynamics simulations. They found clusters diffused by a super-diffusive mechanism (α > 1 in Equation (3)), which is characterized by a continuous displacement of the clusters over the graphene surface without being trapped by adsorption sites. They attributed this behavior to a flat potential energy landscape on the surface, which facilitated the cluster diffusion. Manade et al. [74] estimated the diffusion energy barrier of Cu and Ag adatoms on graphene and arrived at the same conclusion.
Due to their higher adsorption energy, the MSD values of the Cu 6 O 6 cluster are low on both the pristine and functionalized surfaces. Once it was attracted by the surface, it remained stuck to it (see Video S2).   The MSD values corroborate the trends followed by the adsorption energies: the more dynamic the cluster is, the lower the adsorption energy. In agreement with our results, Gervilla et al. [79,80] studied the diffusion of Cu and Ag adatoms and clusters on graphene, using ab initio and classical molecular dynamics simulations. They found clusters diffused by a super-diffusive mechanism (α > 1 in Equation (3)), which is characterized by a continuous displacement of the clusters over the graphene surface without being trapped by adsorption sites. They attributed this behavior to a flat potential energy landscape on the surface, which facilitated the cluster diffusion. Manade et al. [74] estimated the diffusion energy barrier of Cu and Ag adatoms on graphene and arrived at the same conclusion.
Due to their higher adsorption energy, the MSD values of the Cu6O6 cluster are low on both the pristine and functionalized surfaces. Once it was attracted by the surface, it remained stuck to it (see Video S2).

Conclusions
MD simulations were carried out to study the interaction between Cu13, Ag13, and Cu6O6 clusters on pristine and functionalized graphene-based materials with the goal of finding the best structures for adsorption. Adsorption energies, cluster-surface distances, and MSD values show that both the chemical nature of the cluster and the PEGylation of

Conclusions
MD simulations were carried out to study the interaction between Cu 13 , Ag 13 , and Cu 6 O 6 clusters on pristine and functionalized graphene-based materials with the goal of finding the best structures for adsorption. Adsorption energies, cluster-surface distances, and MSD values show that both the chemical nature of the cluster and the PEGylation of the surface are critical to strengthen the interaction between the clusters and surfaces. The introduction of PEG chains favors the adsorption of the three clusters by acting as trapping sites, thus reducing the distances between the clusters and graphene surfaces and increasing the adsorption energy. MSD values also point in this direction. Pristine materials weakly attract Cu 13 and Ag 13 clusters, which keep constantly moving over them. PEG chains act by blocking the displacement of the nanoparticles, which are attracted by their functional groups, thus reducing MSD values. Cu 6 O 6 is more strongly attracted than the metal clusters by both pristine and functionalized materials, as can be inferred from adsorption energies and MSD values. The forcefield used in this work does not allow reproducing the formation of chemical bonds between atoms. Previous ab initio results on metallic clusters deposited on pristine graphene and graphene oxide surfaces show that there may be oxidation of the particles and the formation of covalent bonds, which strengthens the interaction and favors the adsorption process. To correctly describe the behavior of these systems, more accurate quantum mechanical calculations are needed. In addition, the use of an experimentally associated medium (e.g., water) would enable to address in more detail the equilibrium configurations of the structures, as well as the energetic affinity between the clusters and the graphene-based surfaces. Although these more demanding simulations allows gain a deeper insight into the mechanisms that control the interaction between graphene-based materials and nanoparticles, we think our MD results can provide useful qualitative information and allow handling a larger number of atoms in shorter computational times (typical CPU running times of 10-13 h using an Intel ® Xeon ® Quad-core (4 Core) were needed for systems containing 3389 atoms). In future work, DFT calculations will be done to validate these results and be able to consider our molecular dynamics simulations as a quicker reliable method to handle these kinds of systems. Another important aspect that should be considered is the use of a solvent in the calculations. The inclusion of water may have severe effects to several important parameters, such as the equilibrium configurations of the surfaces, the spatial arrangement and the dynamics of the metallic clusters, as well as the energetic affinity between the clusters and the surfaces. Thus, future simulations must include solvent to get a more realistic analysis of real systems.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/nano11061378/s1, Video S1: Cu diffusion on a pristine graphene surface, Video S2: CuO diffusion on a pristine nanoplatelet surface, Video S3: Ag diffusion on a PEGylated nanoplatelet surface, Video S4: Cu cluster movement between two lateral PEG chains in a PEGylated graphene surface.

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article.