Ternary pentagonal BXN (X = C, Si, Ge, and Sn) sheets with high piezoelectricity

Abstract

The discovery of new and stable two-dimensional pentagonal materials with piezoelectric properties is essential for technological advancement. Inspired by recently reported piezoelectric materials penta-BCN and penta-BSiN, we proposed penta-BGeN and penta-BSnN as new members of the penta-family based on first-principles calculations. Comprehensive analyses indicated that both penta-BGeN and penta-BSnN are thermodynamically, dynamically, mechanically, and thermally stable. In terms of mechanical stability, the elastic constant decreased as lower elements in group 4A of the periodic table were used. Therefore, penta-BGeN and penta-BSnN are softer than penta-BCN and penta-BSiN. In terms of piezoelectric properties, piezoelectric stress and strain tensors increase following the same pattern. In group 4A, penta-BSnN had the highest intrinsic piezoelectricity, especially the e₂₂ piezoelectric stress. Typically, the piezoelectric strain d_ij coefficient increases with material softness; penta-BSnN possessed the highest d_ij. Thus, due to its inherent piezoelectricity, penta-BSnN has tremendous potential as a nanoscale piezoelectric material.

---

1 Introduction

Since the discovery of hexagonal graphene, the theoretical exploration and production of two-dimensional (2D) materials have advanced rapidly, even beyond the conventional honeycomb-like structure. A new pentagonal structure, penta-graphene, that resembles the Cairo pentagonal tiling proposed by Zhang et al.,¹ has been extensively investigated through first-principles calculations, although it cannot yet be experimentally synthesized. Its wide band gap of 3.25 eV and high strength could lead to potential applications in nanoelectronics and nanomechanical devices. Other prominent examples are SiC₂, CN₂, BN₂, B₂C, PdS₂, AlN₂, PtN₂, PdN₂, and NiS₂.^2–9

In addition to binary compounds, ternary compounds offer a higher degree of freedom for tuning their electronic and structural properties. The most prominent examples are penta-BCN and hydrogenated penta-BCN,^10,11 which possesses a strong piezoelectric response and is thermodynamically, mechanically, and thermally stable according to density functional theory (DFT). Two dimensional piezoelectric materials have the ability to convert electrical energy into mechanical energy and can find utility in a variety of applications, including sensors, flexible electronics, energy harvesting, and energy storage.^12,13 Subsequently, penta-PdPSe and penta-PdPS were experimentally synthesized.¹⁴ However, penta-PdPSe is non-piezoelectric materials. However, in recent research, the incorporating sulfur into its surface can induce piezoelectricity from non-piezoelectricity in penta-PdPSe.¹⁵

Here, we deployed first-principles calculations to search for a novel member of the ternary pentagonal sheet family with piezoelectric properties. Inspired by the piezoelectricity of penta-BCN, we proposed penta-BGeN and penta-BSnN as two new members of the penta-2D sheet family. We found that both are thermodynamically, dynamically, mechanically, and thermally stable. Penta-BSnN manifests the high piezoelectricity because of its noncentrosymmetric and semiconducting characteristics. Therefore, penta-BSnN is a promising candidate for future nanoscale piezoelectric materials, especially for electromechanical devices.

2 Computational method

First-principles calculations using the framework of DFT are employed with the generalized gradient approximation (GGA) by Perdew–Burke–Ernzerhof (PBE)¹⁶ with the plane-wave basis projector augmented wave (PAW)¹⁷ method, as defined in the Vienna Ab initio Simulation Package (VASP).^18,19 A kinetic energy cutoff of 500 eV is used for the plane-wave expansion. The structures are optimized until the force on each atom is less than 0.01 eV Å⁻¹, and the energy convergence is set to 10⁻⁵ eV. The Monkhorst–Pack²⁰ scheme with 9 × 9 × 1 is employed for the first Brillouin zone integration. A vacuum layer of 20 Å is added along the direction perpendicular to the layer to prevent interactions between a crystal layer and its neighboring image. To verify the dynamic stability, the PHONOPY package²¹ with an energy cutoff of 500 eV and a single k-point at Γ are used to calculate the phonon spectra over a 5 × 5 × 1 supercell. To confirm the thermodynamic stability, ab initio molecular dynamics (AIMD) simulations were deployed under the constant volume and temperature (NVT) ensemble, where the temperature (T) was controlled using the Nose–Hoover thermostat.^22–25 The Heyd–Scuseria–Ernzerhof (HSE06)^26,27 hybrid functional is deployed to calculate the electronic band structure because the standard approximation severely underestimates the band gaps.

3 Results and discussion

We first constructed penta-BXN (X = C, Si, Ge, and Sn) based on the reported atomic configuration of BCN in ref. 10 by replacing the C atoms with Si, Ge, or Sn atoms, as shown in Fig. 1. The pentagonal primitive unit cell, represented by the square in Fig. 1, consists of 6 atoms, including 2 B atoms, 2 X (X = C, Si, Ge, or Sn) atoms, and 2 N atoms. The calculated lattice constants (a and b) and buckling height of penta-BXN (X = C, Si, Ge, and Sn) are listed in Table 1. Our calculated lattice constants for penta-BCN and penta-BSiN matched previous studies.^10,28 The lattice constants a and b increased as we go down a group in the periodic table because the atomic radii increase with elements further down a group. The buckling height h also increases the further down a group an element is.

Fig. 1 3 × 3 supercell of pentagonal BXN (X = C, Si, Ge and Sn) geometrical structures.

Table 1 The calculated lattice constants a and b in Å, bucking height in Å, and the band gaps E^PBE_g and E^HSE_g in eV

	a	b	h	E^PBEg	E^HSEg
penta-BCN	3.670	3.631	1.314	1.70	2.91
penta-BSiN	4.438	4.401	1.602	1.48	2.33
penta-BGeN	4.607	4.590	1.640	1.59	2.42
penta-BSnN	5.045	5.059	1.700	1.05	1.73

---

In binding stability, the cohesive energy per atom was calculated, given by the following equation;^29,30

where E_tot(BXN) is the total energy of the penta-BXN (X = C, Si, Ge and Sn); n_i is the number of Ni, P, and S atoms in the unit cell; and E_i represents the total energies of an isolated single B, X (X = C, Si, Ge or Sn), and N atoms. Our calculated cohesive energies of the penta-BXN (X = C, Si, Ge and Sn) are −7.69, −6.83, −6.24 and −5.84 eV per atom for penta-BCN, penta-BSiN, penta-BGeN and penta-BSnN, respectively. The negative values indicate that all penta-BXN are more energetically favorable than the isolated atoms. We can see that penta-BCN has the lowest cohesive energy per atom. It means that penta-BCN has a strongest binding interaction among penta-BXN. In comparison, the previous report cohesive energy of penta-graphene is −6.84 eV per atom.³¹ For the binary penta-carbides XC₂, the cohesive energy of penta-SiC₂, penta-GeC₂, and penta-SnC₂ are −5.29, −4.69, and −3.83 eV per atom, respectively.³²

The dynamic stability of penta-BCN and penta-BSiN has previously been verified by calculating their phonon dispersion.^10,28 Herein, we verify the dynamic stability of penta-BGeN and penta-BSnN by constructing phonon band structure. The dynamic stability manifests in phonon spectra containing positive frequencies throughout the Brillouin zone. In contrast, if the spectra display imaginary frequency modes (negative values), the structures contain non-restorative forces against a displacement of atoms, indicating dynamic instability. Fig. 2(a) and (b) show the calculated phonon dispersions of penta-BGeN and penta-BSnN, respectively. Thus, as shown in Fig. 2(a) and (b), there are no imaginary modes in the Brillouin zone. These spectra confirm that penta-BXN (X = C, Si, Ge, and Sn) are dynamically stable. We note that the imperceptible negative frequencies at Γ, which are only visible when the spectra are magnified, are a well-known computational error that can be effectively eliminated using the local-density approximations (LDA) functional. These negative frequencies of the phonon dispersion have also been observed in penta-PdSe₂.³³

Fig. 2 Calculated phonon band structures of (a) penta-BGeN and (b) penta-BSnN.

Mechanical stability is an important property for strain engineering. Such stability can be verified from the linear elastic constants, C₁₁, C₂₂, C₁₂, and C₆₆, which are elements in the stiffness tensor. Here, the elastic constants are obtained using the finite difference method. The calculated elastic constants of penta-BXN (X = C, Si, Ge, and Sn) are listed in Table 2. Our calculated elastic constants were comparable to those previously calculated for penta-BCN and penta-BSiN,^10,28 which are listed in Table 2. Typically, structures are considered mechanically stable when the linear elastic constants satisfy the conditions C₁₁C₂₂ – C₁₂² > 0 and C₆₆ > 0. Notably, these conditions were met for all pentagonal BXN, confirming the mechanical stability of penta-BXN (X = C, Si, Ge, and Sn). In comparison, the elastic constant C_ij decreased with increasing atomic radii (i.e., descending a group in the periodic table).

Table 2 The elastic constants C_jk (Nm⁻¹), Poisson's ratio in [100] (ν_a) and [010] (ν_b) directions and Young's modulus (Nm⁻¹) in [100] (E_a) and [010] (E_b) directions

	Previous work		This work
	BCN^10	BSiN^28	BCN	BSiN	BGeN	BSnN
C11	223.56	114.46	222.31	116.97	100.59	66.82
C12	4.90	12.76	4.03	13.06	10.69	17.85
C22	187.39	111.21	187.39	112.60	98.45	66.64
C66	105.38	48.97	105.38	49.36	47.09	32.57
νa	0.022	0.11	0.018	0.112	0.106	0.267
νb	0.026	0.11	0.022	0.116	0.109	0.268
Ea	223.45	113	222.24	115.51	99.46	62.05
Eb	189.03	109	187.30	111.08	97.30	61.86

---

Next, we evaluated the 2D Young's modulus in the [100] and [010] directions (i.e., along the plane), which is obtained by E_a =(C₁₁² − C₁₂²)/C₁₁ and E_b = (C₂₂² − C₁₂²)/C₂₂. Poisson's ratios in the [100] and [010] directions are ν_a = C₁₂/C₁₁ and ν_b = C₁₂/C₂₂, respectively. Young's modulus and Poisson's ratio for all pentagonal BXNs are listed in Table 2. The values of penta-BCN and penta-BSiN are close to those previously reported, E_a = 223.45 Nm⁻¹ and E_b = 189.03 Nm⁻¹ for penta-BCN, and E_a = 113 Nm⁻¹ and E_b = 109 Nm⁻¹ for penta-BSiN. The Young's modulus decreased when we went down a group in the periodic table, while Poisson's ratio increased. These indicate that the pentagonal materials softened as larger elements in group 4A were chosen, and the materials were likely to deform in directions perpendicular to [100] and [010]. Similarly, the elastic constants and Young's modulus of penta-XC₂ (X = C, Si, Ge, and Sn) have been reported.³² The calculated Young's modulus values for penta-XC₂ (X = C, Si, Ge, and Sn) are 266.22, 138.54, 121.65, and 81.23 Nm⁻¹, respectively. It can be observed that in penta-XC₂ (X = C, Si, Ge, and Sn), the Young's modulus decreases as we move down the group in the periodic table. This behavior is similar to that of penta-BXN, where larger elements tend to soften the material. The near-zero Poisson's ratios of penta-BCN indicate that penta-BCN maintains its dimensions under compression and extension along the in-plane directions.

We used eqn (1) and (2) to calculate Young's modulus (E) and Poisson's ratio (ν) along an arbitrary angle θ,^34–36 which is the angle relative to the x direction, in order to further investigate how the mechanical properties of penta-BXN depend on its crystal orientation.

where c = cos θ and s = sin θ. The calculated results for penta-BGeN and penta-BSnN are plotted in Fig. 3, while the calculated results for penta-BCN and penta-BSiN have been provided previously.^10,28 The maxima of the in-plane Young's modulus in penta-BGeN and penta-BSnN are located around the [110] direction (40 for penta-BGeN and 45 for penta-BSnN). Meanwhile, the minima in-plane Young's modulus of both are along [010] (90). In penta-BGeN, the minimum and maximum Poisson's ratios are found along [110] (43) and [010] (90), respectively. In penta-BSnN, the minimum and maximum Poisson's ratios are located at the same angle as penta-BGeN but have higher positive values. Because the in-plane Young's modulus and Poisson's ratios depend on crystal orientation, penta-BXN (X = C, Si, Ge, and Sn) exhibit anisotropic mechanical properties that can be explained by their lattice constants, where a ≠ b.

Fig. 3 Variation of (a) in-plane Young's modulus, and (b) Poisson's ratio for penta-BGeN and penta-BSnN. Note that solid blue and red lines represent penta-BGeN and penta-BSnN, respectively.

Using mapping of indices, the relationship among the elastic (C_jk), piezoelectric stress (e_ik), and strain tensors (d_ij) is

e_ik = d_ijC_jk (4)

Specifically, for a pentagonal system, the relationship between the three variables becomes

Piezoelectric tensors e_ik are calculated using density functional perturbation theory (DFPT) as implemented in the VASP package. The elastic stiffness tensors are calculated by using the finite difference method mentioned before. Therefore, the elements of the strain tensors, d₁₆, d₂₁, d₂₂ and d₃₆, are derived:

According to eqn (5)–(8), the piezoelectric strain tensors d_ij, which are listed in Table 3, are derived from C_jk, as listed in Table 2, and e_ik, as listed in Table 3. Our calculated piezoelectric strain tensors, d_ij, and piezoelectric stress tensors, e_ik, were comparable to those of previous work in penta-BCN and penta-BSiN,^10,28 which are also listed in Table 3. Note that the piezoelectric stress tensors e₂₁ and e₂₂ refer to the induced polarization along the x- and y-directions, respectively, when the strain is applied along the y-direction. Meanwhile, e₁₆ and e₃₆ represent the induced polarization along the xy-direction when the strain is applied along the x- and y-directions, respectively. The negative values of e_ik and d_ij refer to the opposite direction to the strain displacement. From the calculation results, the piezoelectric stress e_ik and strain d_ij tensors tended to increase as the substituted atoms went down the group in the periodic table, which is illustrated in ESI Fig. S1.† It can be inferred that the softer materials possess larger d_ij coefficients. The piezoelectric strain tensor d_ij of penta-BSnN is largest because penta-BSnN possessed the smallest elastic constants. Typically, the piezoelectric strain d_ij coefficient increases with material softness.³⁷ The piezoelectric stress tensor e₂₂ of penta-BSnN is about 3.3 times higher than that of penta-BCN. In addition, all of the piezoelectric coefficients (e₂₁, e₁₆) of penta-BXN are higher than those of two dimensional h-BN (e₁₁ = 1.38 × 10⁻¹⁰ C m⁻¹).³⁸

Table 3 Piezoelectric stress tensors e_ik in unit of × 10⁻¹⁰ C m⁻¹ and piezoelectric strain tensors d_ij of penta-BXN (X = C, Si, Ge and Sn) in the unit of pm V⁻¹

	Previous work		This work
	BCN^10	BSiN^28	BCN	BSiN	BGeN	BSnN
e21	1.93	3.78	1.93	2.35	2.14	1.93
e22	−1.24	−2.40	−1.24	−2.27	−2.68	−4.11
e16	1.80	0.84	1.80	1.50	1.86	2.22
e36	−0.08	−0.26	−0.081	−0.068	−0.113	−0.121
d21	0.878	3.59	0.882	2.27	2.45	4.88
d22	−0.678	−2.57	−0.680	−2.28	−2.98	−7.47
d16	1.72	1.71	1.70	3.03	3.96	6.81
d36	−0.076	−0.53	−0.076	−0.139	−0.240	−0.372

---

The thermal stabilities of penta-BGeN and penta-BSnN were examined via ab initio molecular dynamic simulations. Supercells with a 3 × 3 × 1 duplication of the unit cell are used. The molecular dynamics simulations are performed for 8 ps with a time step of 1 fs at 1200 K. Fig. 4(a) shows the fluctuation of the potential energy of penta-BGeN throughout the simulation at 1200 K. Penta-BGeN is thermally stable at 1200 K. However, penta-BSnN is unable to stand at 1200 K, as shown in ESI Fig. S2.† Therefore, we lower the simulation temperature to 600 K. Finally, penta-BSnN possesses thermal stability around 600 K, as illustrated in Fig. 4(b).

Fig. 4 Total energy fluctuation with time during the AIMD simulation at (a) 1200 K of penta-BGeN and (b) 600 K of penta-BSnN.

The electronic properties of penta-BGeN and penta-BSnN were first explained by the calculated Heyd–Scuseria–Ernzerhof (HSE06)²⁷ band structure, as shown in Fig. 5 When the PBE functional¹⁶ was used, the electronic band gap showed indirect band gaps of 1.59 and 1.05 eV for penta-BGeN and penta-BSnN, as listed in Table 1. The standard DFT underestimates the band gap.^39–41 Therefore, the HSE06 functional,^26,27 which mixes 75% HF exchange and 25% PBE exchange, was also employed to address this issue. With the use of this functional, a larger HSE band gap of 2.42 eV was obtained for penta-BGeN while a larger HSE band gap of 1.73 eV was achieved for penta-BSnN. Both had an indirect band gap character. As shown in ESI Fig. S1,† the valence band maximum (VBM) of penta-BGeN and penta-BSnN are found at Γ point when HSE06 was used, while the conduction band minimum (CBM) was established along the Γ − X path for penta-BGeN, and the along Γ − Y path for penta-BSnN.

Fig. 5 The electronic band structure and partial Density of States (pDOS) of (a) penta-BGeN and (b) penta-BSnN.

Next, we will discuss the effect of band edge position on photocatalytic water splitting. Recent density functional theory (DFT) studies have shown that penta-materials are among the possible candidates for photocatalytic water splitting.^28,32 In water splitting, the band edges position of VBM must cover the redox/oxidation level potentials of water. The standard redox potentials with pH of oxygen evolution reaction (OER) is E_O2/H2O = −5.67 + pH × 0.059 eV, while standard redox potentials with pH of hydrogen evolution reaction (HER) can be express as E_H⁺/H2 = −4.44 + pH × 0.059 eV.^42,43 In the present study, we present the calculated positions of the VBM and CBM with respect to the vacuum level, as shown in ESI Fig. S3.† Our results reveal that the VBM shifts rapidly with increasing pH, while the CBM remains relatively constant. Although, at a pH of 7 (as indicated by the red dashed line), the CBM band edges of all the BXN (X = C, Si, Ge, and Sn) are above the hydrogen evolution potential (H⁺/H₂), only the penta-BCN material has the VBM band edges below the oxygen evolution potential (O₂/H₂O). These findings suggest that penta-BCN is a promising material for water-splitting processes, as it can catalyze both the hydrogen and oxygen evolution reactions at pH 7. Additionally, a decrease in pH (as indicated by the black dashed line) leads to a downward shift in the oxidation and reduction potentials of water, which favors the HER but is detrimental to the OER. Therefore, penta-BCN is the only BXN material suitable for OER under highly acidic conditions (pH = 0) in water-splitting processes.

The band gaps of 2D materials can potentially be modulated by strain. Thus, we also investigated how the band gaps of penta-BGeN and penta-BSnN change with additional biaxial strain from equilibrium (i.e., zero strain) to the critical strain point. In order to avoid the interaction of periodic conditions, the 2 × 2 × 1 supercell has been adopted to release the special constraints. The values of two critical strains, (ε_C1 and ε_C2), can be inferred from Fig. 6. The first is the point where the derivative curve reaches its maximum turning point. This point is widely known as the bond-breaking point. The second is the point at which the plastic range begins to manifest irreversible deformations (reconstruction process). Above ε_C2, the strain energy always increases, and the system maintains its reconstructive structure. In penta-BGeN, ε_C1 occurred at nearly ε = 0.14 while it occurred at nearly ε = 0.08 in penta-BSnN. This means that penta-BGeN could withstand biaxial strain before mechanical failure at these critical points. However, in penta-BGeN, structural instability occurred at ε_C2, at which point the structure starts to collapse.

Fig. 6 Strain–stress relations for (a) penta-BGeN and (b) penta-BSnN with biaxial strains. The vertical dashed lines indicate the critical structure instability.

4 Conclusions

By using DFT-based calculations, we theoretically explored penta-BGeN and penta-BSnN as the new members in the penta-2D family. Both penta-BGeN and penta-BSnN are dynamically, mechanically, and thermally stable, as comprehensively verified by the phonon dispersion, elastic constants, and molecular dynamics simulations, respectively. Penta-BGeN and penta-BSnN are semiconductors with indirect band gaps of 2.42 and 1.73 eV, respectively. Penta-BGeN and penta-BSnN are soft materials because their elastic constants decrease as we go down the group in the periodic table. The elements of the piezoelectric stress and strain tensors increase as we move down a group in the periodic table. In the same family, penta-BSnN yields the highest intrinsic piezoelectricity, especially the e₂₂ piezoelectric stress. Typically, the piezoelectric strain, d_ij, coefficient increases with material softness, yielding the greatest d_ij for penta-BSnN. Thus, penta-BSnN has promising prospects as a nanoscale piezoelectric material due to its inherent piezoelectricity.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This research has received funding support from the NSRF via the Program Management Unit for Human Resources & Institutional Development, Research and Innovation [Grant Number B05F650024]. T. T. was supported by Graduate Program Scholarship from the Graduate School, Kasetsart University. Proofread and any publication fee are supported by International SciKU Branding (ISB), Faculty of Science, Kasetsart University. We wish to thank NSTDA Supercomputer Center (ThaiSC) for providing computing resources for this work.

References

1. S. Zhang, J. Zhou, Q. Wang, X. Chen, Y. Kawazoe and P. Jena, Proc. Natl. Acad. Sci. U. S. A., 2015, 112, 2372–2377.
2. G. R. Berdiyorov and M. E.-A. Madjet, RSC Adv., 2016, 6, 50867–50873.
3. S. Zhang, J. Zhou, Q. Wang and P. Jena, J. Phys. Chem. C, 2016, 120, 3993–3998.
4. J. Li, X. Fan, Y. Wei and G. Chen, Sci. Rep., 2016, 6, 31840.
5. F. Li, K. Tu, H. Zhang and Z. Chen, Phys. Chem. Chem. Phys., 2015, 17, 24151–24156.
6. Z. Cheng, X. Zhang, H. Zhang, J. Gao, H. Liu, X. Yu, X. Dai, G. Liu and G. Chen, Phys. Chem. Chem. Phys., 2021, 23, 6278–6285.
7. Y.-S. Lan, X.-R. Chen, C.-E. Hu, Y. Cheng and Q.-F. Chen, J. Mater. Chem. A, 2019, 7, 11134–11142.
8. J. Li, X. Fan, Y. Wei, H. Liu, S. Li, P. Zhao and G. Chen, Sci. Rep., 2016, 6, 1–10.
9. C.-T. Wang and S. Du, Phys. Chem. Chem. Phys., 2020, 22, 7483–7488.
10. K. Zhao, Y. Guo, Y. Shen, Q. Wang, Y. Kawazoe and P. Jena, J. Phys. Chem. Lett., 2020, 11, 3501–3506.
11. T. Thanasarnsurapong, K. Dabsamut, T. Maluangnont, J. T-Thienprasert, S. Jungthawan and A. Boonchun, J. Appl. Phys., 2021, 129, 095101.
12. S. D. Mahapatra, P. C. Mohapatra, A. I. Aria, G. Christie, Y. K. Mishra, S. Hofmann and V. K. Thakur, Adv. Sci., 2021, 8, 2100864.
13. M. B. Ghasemian, T. Daeneke, Z. Shahrbabaki, J. Yang and K. Kalantar-Zadeh, Nanoscale, 2020, 12, 2875–2901.
14. P. Li, J. Zhang, C. Zhu, W. Shen, C. Hu, W. Fu, L. Yan, L. Zhou, L. Zheng and H. Lei, et al., Adv. Mater., 2021, 33, 2102541.
15. K. Dabsamut, I. Chatratin, T. Thanasarnsurapong, S. Jungthawan and A. Boonchun, J. Alloys Compd., 2023, 169640.
16. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
17. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953.
18. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
19. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169.
20. H. J. Monkhorst and J. D. Pack, Phys. Rev. B: Condens. Matter Mater. Phys., 1976, 13, 5188.
21. A. Togo and I. Tanaka, Scr. Mater., 2015, 108, 1–5.
22. S. Nosé, J. Chem. Phys., 1984, 81, 511–519.
23. N. Shuichi, Prog. Theor. Phys. Suppl., 1991, 103, 1–46.
24. W. G. Hoover, Phys. Rev. A, 1985, 31, 1695.
25. D. Frenkel and B. Smit, Understanding molecular simulation: from algorithms to applications, Elsevier, 2001, vol. 1.
26. J. Heyd, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys., 2003, 118, 8207–8215.
27. J. Heyd, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys., 2006, 124, 219906.
28. M. J. Varjovi, M. Kilic and E. Durgun, Phys. Rev. Mater., 2022, 6, 034004.
29. K. Dabsamut, T. Thanasarnsurapong, I. Chatratin, T. Maluangnont, S. Jungthawan and A. Boonchun, J. Phys. Chem. C, 2022, 126, 19455–19461.
30. K. Dabsamut, I. Chatratin, T. Thanasarnsurapong, T. Maluangnont and A. Boonchun, Phys. Chem. Chem. Phys., 2023, 25, 3815–3819.
31. N. Sathishkumar, S.-Y. Wu and H.-T. Chen, Chem. Eng. J., 2020, 391, 123577.
32. M. E. Kilic and K.-R. Lee, Phys. Rev. Mater., 2021, 5, 065404.
33. A. D. Oyedele, S. Yang, L. Liang, A. A. Puretzky, K. Wang, J. Zhang, P. Yu, P. R. Pudasaini, A. W. Ghosh and Z. Liu, et al., J. Am. Chem. Soc., 2017, 139, 14090–14097.
34. E. Cadelano, P. L. Palla, S. Giordano and L. Colombo, Phys. Rev. B, 2010, 82, 235414.
35. H. Wang, X. Li, P. Li and J. Yang, Nanoscale, 2017, 9, 850–855.
36. Y. Zhao, X. Li, J. Liu, C. Zhang and Q. Wang, J. Phys. Chem. Lett., 2018, 9, 1815–1820.
37. R. Hinchet, U. Khan, C. Falconi and S.-W. Kim, Mater. Today, 2018, 21, 611–630.
38. K.-A. N. Duerloo, M. T. Ong and E. J. Reed, J. Phys. Chem. Lett., 2012, 3, 2871–2876.
39. K. Dabsamut, T. Maluangnont, P. Reunchan, J. T-Thienprasert, S. Jungthawan and A. Boonchun, Appl. Phys. Lett., 2022, 120, 203101.
40. H. Xiao, J. Tahir-Kheli and W. A. Goddard III, J. Phys. Chem. Lett., 2011, 2, 212–217.
41. M. K. Chan and G. Ceder, Phys. Rev. Lett., 2010, 105, 196403.
42. G. W. VanLoon and S. J. Duffy, Environmental chemistry: a global perspective, Oxford university press, 2017.
43. C. Dues, W. G. Schmidt and S. Sanna, ACS Omega, 2019, 4, 3850–3859.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d2ra08342f

---