Stacking interactions of nickel bis(dithiolene) with graphene and beyond

Abstract

Density functional theory (DFT) with dispersion correction has been used to study the stacking interactions of the nickel bis(dithiolene) molecule and graphene. As in our previous study of the nickel bis(dithiolene) molecule and benzene, two different configurations were considered for the nickel bis(dithiolene) molecule and graphene, and the whole potential energy surface (PES) was explored for each of them. The stacking interaction energy is comparable with other similar systems in the experiments, and it is shown that the surface-mediated interactions play a role in determining molecular orientation. Based on the results of the nickel bis(dithiolene) molecule and graphene, we also studied a new two-dimensional (2D) heterobilayered material consisting of a 2D nickel bis(dithiolene) sheet and graphene or hexagonal boron nitride (h-BN). It turns out that h-BN is a good substrate for the 2D nickel bis(dithiolene) sheet, an organic topological insulator, to maintain its original properties, while the 2D nickel bis(dithiolene) sheet behaves like a metal when forming a new 2D heterobilayered material with graphene. These observations would undoubtedly enrich the current study of 2D heterobilayered materials.

---

The single-atom thick graphene has attracted tremendous attention in various disciplines for both fundamental interest^1,2 and many possible applications, such as gas sensors,³ hydrogen storage,⁴ polymer–graphene nanocomposites,⁵ plasmonic devices,⁶ membranes,⁷ and coatings.⁸ Following graphene, other 2D layered materials, such as h-BN, transition metal dichalcogenides (e.g. MoS₂ and NbSe₂), and group V metal chalcogenides (e.g. Sb₂Te₃ and Bi₂Se₃), have since been prepared by exfoliation and chemical vapor deposition.⁹ These new 2D materials revealed a wealth of novel physics including 2D gas of massless Dirac fermions, ambipolar electric field effect, room-temperature quantum Hall effect, breakdown of the Born–Oppenheimer approximation, and quantum capacitance.^1,10–13 However, it is until recently that great efforts have been made for tuning the properties of the layered materials by creating hybrid structures. As one kind of hybrid structures, the heterobilayered 2D materials are usually van der Waals (vdW) bonded, and show emergent properties and functionalities of the individual layers.^9,14 Due to the excellent lattice match of graphene and h-BN, the heterobilayered graphene/h-BN 2D materials have at first been investigated both theoretically and experimentally.^15–19 The most stable configuration for heterobilayered graphene/h-BN is that one C atom lies on top of a B atom, while the other C atom sits above the center of the h-BN ring. Besides the heterobilayered graphene/h-BN 2D materials, there are also a number of multi-layered structures that have been proposed theoretically. For instance, the bilayers of MoS₂/graphene and MoS₂/WS₂ have been demonstrated by theory to have a potential for solar energy absorption and conversion at the nanoscale.²⁰

To date, most of the 2D layered materials are generated from exfoliation of a 3D layered solid, which is a top-down method. Nevertheless, Nishihara et al. recently have synthesized a highly π-conjugated 2D nickel bis(dithiolene) sheet by a bottom-up method using benzenehexathiol (BHT) and nickel(II) acetate (Ni(OAc)₂).²¹ This 2D organometallic framework, adopting hexagonal symmetry, is suggested to be an organic topological insulator (OTI) by subsequent calculations.²² In addition, it is also found that the ethylene molecules can be bound and released by the 2D nickel bis(dithiolene) sheet under different conditions,²³ which is similar to nickel bis(dithiolene) molecule. Nickel bis(dithiolene) molecule is a chelate, and has been widely studied due to its unusual properties.^24–26 Recent experiments demonstrated that graphene has a strong interaction with copper phthalocyanine (CuPc), and could control the orientation of CuPc on its surface.²⁷ The highly ordered organic crystals on graphene have a great potential in the device for charge transport, exciton diffusion, and dissociation.²⁸

Our previous calculations have shown that nickel bis(dithiolene) molecule has strong interaction with benzene molecule.²⁹ In the current work, we will use DFT methods to study the interactions of the nickel bis(dithiolene) molecule on graphene, searching for the most favorable orientation. Based on the results of the nickel bis(dithiolene) molecule on graphene and also inspired by the recent developments of 2D heterobilayered materials, we continue to investigate the feasibility of 2D nickel bis(dithiolene)/graphene and nickel bis(dithiolene)/h-BN heterobilayers, as well as their electronic properties.

First-principles calculations were carried out using the Vienna ab initio simulation package (VASP).^30–33 The Kohn–Sham equations were solved using the projector-augmented wave (PAW) method.^34,35 The exchange and correlation interactions of valence electrons are described by the Perdew–Burke–Ernzerhof (PBE) functional^36,37 within the generalized gradient approximation (GGA). Since vdW interactions are expected to be significant in the stacking systems, DFT-D2 method of Grimme³⁸ was also applied throughout. The Brillouin-zone integrations were performed on a Γ-centered Monkhorst–Pack 4 × 4 × 1 k-point grid³⁹ for 2D structures, and Γ point only for the nickel bis(dithiolene) molecule and 30 × 30 Å supercell graphene. The kinetic energy cutoff for plane waves was set to 400 eV and the “accurate” precision setting was adopted to avoid wrap around errors. The convergence criterion for the electronic self-consistent loop was set to 10⁻⁵ eV. During the structure optimizations, the vacuum regions were at least 20 Å to ensure the periodic images are well separated while other lattice vectors were fully relaxed. All atoms were also relaxed until the Hellmann–Feynman forces were smaller than 0.01 eV Å⁻¹.

First, we studied the stacking interaction of the nickel bis(dithiolene) molecule, Ni(S₂C₆H₄)₂, shown in Fig. 1(a) on top of graphene. Both the nickel bis(dithiolene) molecule and graphene (30 × 30 Å supercell is used to ensure that the Ni complexes on the top are well separated with ca. 22.3 Å between each other) are optimized separately. The optimized Ni–S bond length is 2.137 Å and S–C bond length is 1.717 Å. These values are close to the counterparts of the nickel bis(dithiolene) framework of our study (Ni–S is 2.135 Å and C–S is 1.706 Å), seen in Fig. 1(b), and of the others (Ni–S is 2.138 Å and C–S is 1.707 Å),²³ indicating that the molecular nickel bis(dithiolene) well resembles its 2D solid form. During the calculation of the stacking interaction, only the z components (graphene is on XY plane) of the atoms of nickel bis(dithiolene) molecule are allowed to move, while all the others are kept fixed.

Fig. 1 The optimized structures of the nickel bis(dithiolene) molecule (a) and the 2D nickel bis(dithiolene) sheet (b). The solid rhombus demonstrates the unit cell of the 2D sheet. The colors of the atoms are kept the same thereafter.

Similar to our previous work,²⁹ two orientations of the nickel bis(dithiolene) molecule with respect to graphene have been considered in this report, and they are different by 90 degree rotation, as shown in Fig. 2(a) and (c) for the respective top views. The side views of configurations A and B are similar, and Fig. 2(e) represents both of them with the normal distance H. In order to search the most stable position of the nickel bis(dithiolene) molecule on top of graphene, we explored the PES for these two configurations, as the nickel bis(dithiolene) molecule moves over graphene. The red trapezoids in Fig. 2(b) and (d) demonstrate the minimum repeated areas in one six-membered ring of graphene for A and B configurations, respectively. The displacement of the Ni atom of the nickel bis(dithiolene) molecule relative to graphene on X and Y axes is measured within the red trapezoidal areas and interaction energy E is calculated by formula (1), with the corresponding normal distance H. Based upon the symmetry of graphene and the nickel bis(dithiolene) molecule, the whole PES can be derived accordingly.

E = E(nickel bis(dithiolene)/graphene complex) − E(nickel bis(dithiolene)) − E(graphene) (1)

Fig. 2 The molecular nickel bis(dithiolene)/graphene stacking systems used for calculations: top view (a) for configuration A and top view (c) for configuration B. The displacement of the Ni atom of the molecular nickel bis(dithiolene) relative to graphene on X and Y axes is measured within the red trapezoidal areas in (b) and (d) for A and B configurations, respectively. The x label marks the most stable position of the Ni atom in each case. (e) shows the side view with the normal distance (H) between the nickel bis(dithiolene) molecule and graphene.

Fig. 3 shows the calculated normal distance H and interaction energy E as a function of X and Y for the molecular nickel bis(dithiolene)/graphene configuration A in one six-membered ring of graphene. After optimization, the nickel bis(dithiolene) molecule is no longer on the same plane, but tilts a little toward graphene. Thus, the optimal distance H is set to the distance of the Ni atom relative to graphene. The optimal distance H reaches the global minimum, 3.41 Å, when the Ni atom sitting right on the C atoms of graphene, and the global maximum, 3.47 Å, when the Ni atom sitting on the center of six-membered ring of graphene. The minimum value of the optimal H is close the counterpart (ca. 3.3 Å at the PBE-vdW level) of CuPc on graphene in recent experiments,²⁷ and both of them have a metal–carbon head to head configuration. The interaction energy also reaches the global maximum (−1.39 eV) when the optimal distance H is the maximum, seen in Fig. 3(c) and (d), but the global minimum (−1.48 eV) happens when the Ni atom is right on top of the middle of certain C–C bonds of graphene (H = 3.42 Å), as is labeled as x in Fig. 2(b). The interaction energy is ca. one half of that of CuPc on graphene (−3.37 eV at the PBE-vdW level),²⁷ even though CuPc is much bigger in size than the nickel bis(dithiolene) molecule. In addition, it should also be pointed out that the interaction energy cannot be directly compared due to the different treatments of the dispersion correction. Similarly, Fig. 4 shows the calculated normal distance H and interaction energy E as a function of X and Y for the molecular nickel bis(dithiolene)/graphene configuration B. However, it is interesting that for configuration B, the interaction energy reaches the global minimum (−1.44 eV) when the optimal distance H is the maximum (3.48 Å), and the Ni position is labeled as x in Fig. 2(d) as well. Fig. 2(c) shows the top view of the most stable configuration B. The most stable configuration B is slightly higher in energy than its counterpart of configuration A, by 0.04 eV (less than 1.0 kcal mol⁻¹), which is similar to our previous study on the molecular nickel bis(dithiolene)/benzene.²⁹ Our calculations mimic a low coverage of nickel bis(dithiolene) molecules on graphene. It could be concluded that under the low coverage condition, nickel bis(dithiolene) molecule prefers configuration A on graphene. As CuPc on graphene,²⁷ one could expect to see surface-induced orientation control of nickel bis(dithiolene) molecules on graphene as well.

Fig. 3 The calculated normal distance H (in Å) (3D view (a) and top view (b)) and interaction energy E (in eV) (3D view (c) and top view (d)) as a function of X and Y for the molecular nickel bis(dithiolene)/graphene configuration A. The origin is set to the center of one six-membered ring of graphene.

Fig. 4 The calculated normal distance H (in Å) (3D view (a) and top view (b)) and interaction energy E (in eV) (3D view (c) and top view (d)) as a function of X and Y for the molecular nickel bis(dithiolene)/graphene configuration B. The origin is set to the center of one six-membered ring of graphene.

Now we turn our attention to the 2D nickel bis(dithiolene)/graphene heterobilayer. The geometry of the 2D nickel bis(dithiolene) sheet is shown Fig. 1(b), along with the unit cell, and its band structure is shown in Fig. 6(c). Both the geometric and electronic properties are in excellent agreement with the previous studies.^22,23 The 2D nickel bis(dithiolene) sheet has a same hexagonal symmetry with graphene. The optimized lattice parameter of the 2D nickel bis(dithiolene) sheet unit cell is 14.60 Å, while the lattice parameter of graphene unit cell is 2.46 Å at the same level. The unit cell of 2D nickel bis(dithiolene) is almost six times that of graphene, with a mismatch ca. 1.0%. This mismatch is less than that of graphene with h-BN (2%),¹⁵ and of graphene with MoS₂ (1.3%).²⁰ Thus, a unit 2D nickel bis(dithiolene) was placed on top of a 6 × 6 hexagonal supercell of graphene to form a heterobilayer, as shown in Fig. 5(a)–(c). This configuration comes from the most stable configuration B, because configuration A does not match the 2D nickel bis(dithiolene) in symmetry.

Fig. 5 The optimized unit cell of the most stable 2D nickel bis(dithiolene)/graphene sheet: top view by spacefill model (a); side view by spacefill model (b) with normal distance R. The top view by ellipsoid model for the 2D nickel bis(dithiolene)/graphene (c) and nickel bis(dithiolene)/h-BN sheet (d). (e) demonstrates the calculated interaction energy (in eV) as a function of R for the 2D nickel bis(dithiolene)/graphene (Ni_C) and nickel bis(dithiolene)/h-BN (Ni_BN) sheet.

The new 2D heterobilayer adopts a unit cell, with lattice constant 14.74 Å. The interaction energy as a function of the distance R between two layers, as shown in Fig. 5(b), is plotted in Fig. 5(e). The energy at the minimum is −2.02 eV, and the corresponding R is 3.45 Å. The calculated minimum interaction energy is much lower than that of graphene/h-BN counterpart (6 × 6 supercell), while the layer–layer distance is almost the same.¹⁷ Recent experiments have made the graphene/h-BN heterobilayer available,^16,18,19 and our proposed 2D heterobilayer could be another good candidate for synthesis, and expand the current graphene-based heterobilayer study. Because of the similarity of graphene and h-BN in geometry, we also studied the 2D nickel bis(dithiolene)/h-BN heterobilayer, as seen in Fig. 5(d). The interaction energy of the 2D nickel bis(dithiolene) and h-BN reaches the minimum, −2.03 eV, when R is 3.35 Å. It appears that the 2D nickel bis(dithiolene)/h-BN heterobilayer is more stable than its nickel bis(dithiolene)/graphene counterpart. The 2D nickel bis(dithiolene) sheet is an organic topological insulator, while graphene is a semi-metal, and h-BN is a good insulator. It is interesting to see how the electronic properties of the 2D nickel bis(dithiolene) sheet change when it is put on the graphene or h-BN. Fig. 6(a) and (b) show the band structures of the 2D nickel bis(dithiolene)/graphene and nickel bis(dithiolene)/h-BN heterobilayers, respectively. For nickel bis(dithiolene)/h-BN heterobilayer, it is clear to see that the valence band and conduction band of it comes exactly from the 2D nickel bis(dithiolene) sheet, which is understandable since h-BN is a good insulator, with a wide band gap (ca. 5 eV). However, when forming new 2D heterobilayered material with graphene, the conduction band of the 2D nickel bis(dithiolene) sheet goes way down under the Fermi level, assuming a metallic behavior as seen in Fig. 6(a). This phenomenon has something to do with the nature of graphene, whose well delocalized π electrons can transfer to the conduction band of the 2D nickel bis(dithiolene) sheet. The dramatic difference of the two heterobilayers demonstrates two unique ways to utilize the newly synthesized 2D nickel bis(dithiolene) sheet: as a good substrate, h-BN could be used to sustain the 2D nickel bis(dithiolene) sheet with a strong interaction but maintain its electronic properties, while graphene could easily tune the properties of the 2D nickel bis(dithiolene) sheet, from insulator to metal, with a strong interaction as well.

Fig. 6 Band structures of the 2D nickel bis(dithiolene)/graphene sheet (a), the 2D nickel bis(dithiolene)/h-BN sheet (b), and the 2D nickel bis(dithiolene) sheet (c).

In summary, we use DFT with dispersion correction to study the stacking interactions of nickel bis(dithiolene) with graphene. We find the most stable configuration of the nickel bis(dithiolene) molecule on graphene by mapping the whole PES. The stacking interaction between the nickel bis(dithiolene) molecule and graphene is very strong, comparable with the previous experiments of CuPc with graphene, indicating that there would be surface-induced orientation for the nickel bis(dithiolene) molecule on graphene. In light of what we have found in the nickel bis(dithiolene) molecule and graphene, we continue to explore the possibility of a new 2D heterobilayered material consisting of the 2D nickel bis(dithiolene) sheet with graphene and h-BN. Given the perfect match of the lattice constant as well as the same crystal symmetry, the 2D nickel bis(dithiolene) sheet and graphene or h-BN would be an ideal model for graphene-based organic electronics, and compensate the current study on the heterobilayer of graphene and h-BN. The results suggest that the 2D nickel bis(dithiolene) sheet and graphene or h-BN can form a stable heterobilayered material, with the layer–layer distance close to that of graphene and h-BN, but the electronic properties of the two heterobilayers are quite different. We are confident that these calculations will stimulate experimental studies on the stacking interaction between nickel bis(dithiolene) and graphene or h-BN.

This work used computational resources of the National Center for Computational Sciences at Oak Ridge National laboratory and of the National Energy Research Scientific Computing Center, which are supported by the Office of Science of the U.S. Department of Energy under Contract no. DE-AC05-00OR22750 and DE-AC02-05CH11231, respectively. We also acknowledge the support from the Center for Nanophase Materials Sciences, which is sponsored at ORNL by the Scientific User Facilities Division, U.S. Department of Energy.

References

1. K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos and A. A. Firsov, Nature, 2005, 438, 197–200.
2. K. S. Novoselov, E. McCann, S. V. Morozov, V. I. Fal'ko, M. I. Katsnelson, U. Zeitler, D. Jiang, F. Schedin and A. K. Geim, Nat. Phys., 2006, 2, 177–180.
3. R. Arsat, M. Breedon, M. Shafiei, P. G. Spizziri, S. Gilje, R. B. Kaner, K. Kalantar-zadeh and W. Wlodarski, Chem. Phys. Lett., 2009, 467, 344–347.
4. A. Politano and G. Chiarello, Carbon, 2013, 61, 412–417.
5. V. Alzari, V. Sanna, S. Biccai, T. Caruso, A. Politano, N. Scaramuzza, M. Sechi, D. Nuvoli, R. Sanna and A. Mariani, Compos. Part B: Eng., 2014, 60, 29–35.
6. A. Politano and G. Chiarello, Nanoscale, 2013, 5, 8215–8220.
7. V. Berry, Carbon, 2013, 62, 1–10.
8. M. Topsakal, H. Şahin and S. Ciraci, Phys. Rev. B, 2012, 85, 155445.
9. S. Z. Butler, S. M. Hollen, L. Cao, Y. Cui, J. A. Gupta, H. R. Gutiérrez, T. F. Heinz, S. S. Hong, J. Huang, A. F. Ismach, E. Johnston-Halperin, M. Kuno, V. V. Plashnitsa, R. D. Robinson, R. S. Ruoff, S. Salahuddin, J. Shan, L. Shi, M. G. Spencer, M. Terrones, W. Windl and J. E. Goldberger, ACS Nano, 2013, 7, 2898–2926.
10. K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva and A. A. Firsov, Science, 2004, 306, 666–669.
11. Y. Zhang, Y.-W. Tan, H. L. Stormer and P. Kim, Nature, 2005, 438, 201–204.
12. S. Pisana, M. Lazzeri, C. Casiraghi, K. S. Novoselov, A. K. Geim, A. C. Ferrari and F. Mauri, Nat. Mater., 2007, 6, 198–201.
13. J. Xia, F. Chen, J. Li and N. Tao, Nat. Nanotechnol., 2009, 4, 505–509.
14. Q. Tang and Z. Zhou, Prog. Mater. Sci., 2013, 58, 1244–1315.
15. G. Giovannetti, P. A. Khomyakov, G. Brocks, P. J. Kelly and J. van den Brink, Phys. Rev. B, 2007, 76, 073103.
16. R. g. Decker, Y. Wang, V. W. Brar, W. Regan, H.-Z. Tsai, Q. Wu, W. Gannett, A. Zettl and M. F. Crommie, Nano Lett., 2011, 11, 2291–2295.
17. Y. Fan, M. Zhao, Z. Wang, X. Zhang and H. Zhang, Appl. Phys. Lett., 2011, 98, 083103.
18. J. Xue, J. Sanchez-Yamagishi, D. Bulmash, P. Jacquod, A. Deshpande, K. Watanabe, T. Taniguchi, P. Jarillo-Herrero and B. J. LeRoy, Nat. Mater., 2011, 10, 282–285.
19. W. Yang, G. Chen, Z. Shi, C.-C. Liu, L. Zhang, G. Xie, M. Cheng, D. Wang, R. Yang, D. Shi, K. Watanabe, T. Taniguchi, Y. Yao, Y. Zhang and G. Zhang, Nat. Mater., 2013, 12, 792–797.
20. M. Bernardi, M. Palummo and J. C. Grossman, Nano Lett., 2013, 13, 3664–3670.
21. T. Kambe, R. Sakamoto, K. Hoshiko, K. Takada, M. Miyachi, J.-H. Ryu, S. Sasaki, J. Kim, K. Nakazato, M. Takata and H. Nishihara, J. Am. Chem. Soc., 2013, 135, 2462–2465.
22. Z. F. Wang, N. Su and F. Liu, Nano Lett., 2013, 13, 2842–2845.
23. Q. Tang and Z. Zhou, J. Phys. Chem. C, 2013, 117, 14125–14129.
24. P. Cassoux, L. Valade, H. Kobayashi, A. Kobayashi, R. A. Clark and A. E. Underhill, Coord. Chem. Rev., 1991, 110, 115–160.
25. R. K. Szilagyi, B. S. Lim, T. Glaser, R. H. Holm, B. Hedman, K. O. Hodgson and E. I. Solomon, J. Am. Chem. Soc., 2003, 125, 9158–9169.
26. P. Deplano, L. Pilia, D. Espa, M. L. Mercuri and A. Serpe, Coord. Chem. Rev., 2010, 254, 1434–1447.
27. K. Xiao, W. Deng, J. K. Keum, M. Yoon, I. V. Vlassiouk, K. W. Clark, A.-P. Li, I. I. Kravchenko, G. Gu, E. A. Payzant, B. G. Sumpter, S. C. Smith, J. F. Browning and D. B. Geohegan, J. Am. Chem. Soc., 2013, 135, 3680–3687.
28. A. Schlierf, P. Samori and V. Palermo, J. Mater. Chem. C, 2014 10.1039/c3tc32153c.
29. J. Zhou, Chem. Phys. Lett., 2013, 591, 29–31.
30. G. Kresse and J. Hafner, Phys. Rev. B, 1993, 47, 558–561.
31. G. Kresse and J. Hafner, Phys. Rev. B, 1994, 49, 14251–14269.
32. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
33. G. Kresse and J. Furthmüller, Phys. Rev. B, 1996, 54, 11169–11186.
34. P. E. Blöchl, Phys. Rev. B, 1994, 50, 17953–17979.
35. G. Kresse and D. Joubert, Phys. Rev. B, 1999, 59, 1758–1775.
36. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
37. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1997, 78, 1396.
38. S. Grimme, J. Comput. Chem., 2006, 27, 1787–1799.
39. H. J. Monkhorst and J. D. Pack, Phys. Rev. B, 1976, 13, 5188–5192.

---