INTRODUCTION
Environmental pollution, especially water pollution, has become one of the most pressing global challenges in the 21st century 1. Industrial growth, agricultural waste, and massive urbanization have significantly increased hazardous pollutants in natural water systems, including heavy metals, dyes, and other organic compounds 2. These pollutants pose serious risks to human health and disrupt the balance of ecosystems 3. Consequently, there is an urgent need to develop efficient and sustainable materials for water purification. Among the various methods explored, adsorption has high effectiveness, operational simplicity, and relatively low cost 4. However, adsorption is highly dependent on the characteristics of the adsorbent material, highlighting the need for materials with high selectivity and stability 5.
Graphene, a two-dimensional carbon, has attracted great attention in the environmental field due to its extraordinary physical and chemical properties 6. Its flat and ultrathin structure, very large specific surface area, high mechanical strength, and excellent electrical conductivity make graphene an ideal candidate for adsorption applications 7. In addition, graphene surfaces can be chemically modified to enhance interactions with polar molecules such as water (H2O) and dissolved contaminants 8. Recent studies have demonstrated the selective adsorption ability of graphene to water molecules, which opens up new possibilities in the development of nanotechnology materials for water treatment 9.
Although the potential of graphene is growing in environmental applications, accurate modeling of the interaction between graphene and adsorbates requires precise computational parameters, especially in density functional theory (DFT) simulations 10. One of the important parameters in DFT calculations is the k-point mesh, which refers to the selection of sampling points in the Brillouin zone of the periodic system 11. Inappropriate or too coarse k-point selection can lead to incorrect predictions of adsorption energy, electronic properties, and structural stability 12. Therefore, this study systematically investigates the effect of k-point mesh density on the adsorption behavior of H2O molecules on graphene surfaces 13. By emphasizing the important role of k-point selection, this study aims to establish an optimal computational framework for future theoretical studies on graphene–adsorbate interactions.
METHODS
We employed first-principles density functional theory (DFT) calculations to study the adsorption behavior of H2O molecules on graphene. The simulations were performed using Quantum Espresso 14, utilizing projector-augmented wave (PAW) pseudopotentials. Previous study already showed that Quantum Espresso gave adsorption energies description of the sodium ion calculation15. The wavefunction cutoff energy was set to 60 Ry, while the charge density cutoff energy was defined at 600 Ry. To determine the best k-point to sampling the Brillouin zone, we used a (8 × 8 × 1) until (20 × 20 × 1) k-point mesh with offset/ No offset mesh. A (1×1) supercell was used in this calculation with 15 Å vacuum slabs. We also used the DFT-D3 functional as an exchange-correlation functional. This functional is used to account for the van der Waals (vdW) interactions, which is used to describe the weak interactions in the system. All the calculations were performed until all the forces acting to each atom in the system were less than 10-3 Ry. We then compared each energy with other k-points mesh variation to obtain the appropriate for our calculations.
| k-Points | Offset/No Offset | Lattice A | Lattice B | C–C Bond Length (Å) | Calculated Total Energy (Ry) | Computational Time |
|---|---|---|---|---|---|---|
| 2×8×1 | Offset | 2.4651 | 2.4643 | 1.4231 | −24.8628 | 11m 49.18s |
| No Offset | 2.4651 | 2.4648 | 1.4234 | −24.8627 | 9m 2.32s | |
| 2×12×1 | Offset | 2.4613 | 2.4665 | 1.4206 | −24.8623 | 5m 11.86s |
| No Offset | 2.4650 | 2.4645 | 1.4232 | −24.8628 | 7m 48.31s | |
| 2×16×1 | Offset | 2.4646 | 2.4647 | 1.4240 | −24.8628 | 20m 4.30s |
| No Offset | 2.4646 | 2.4646 | 1.4260 | −24.8628 | 10m 46.75s | |
| 2×20×1 | Offset | 2.4648 | 2.4648 | 1.4233 | −24.8627 | 24m 3.95s |
| No Offset | 2.4657 | 2.4647 | 1.4232 | −24.8627 | 8m 23.42s |
In order to obtain the total energy, we employ the Kohn–Sham formula 16 of Density Functional Theory, where the interacting many-electron problem is replaced by a system of non-interacting electrons governed by
Here, Ts[ρ] represents the kinetic energy of the non-interacting electrons, Eext[ρ] corresponds to the interaction energy between electrons and the external potential (typically arising from nuclei), EH[ρ] is the classical electrostatic (Hartree) energy, and Exc[ρ] is the exchange-correlation energy, which includes all many-body interactions beyond the Hartree approximation.
The Kohn–Sham equations are solved self-consistently to obtain the ground-state electron density ρ(r). In our calculations, exchange-correlation effects are treated using the Perdew–Burke–Ernzerhof (PBE) functional within the generalized gradient approximation (GGA). PAW pseudopotentials represent core electrons, and valence electrons are expanded in a plane-wave basis set with a cutoff energy of 60 Ry. Structural optimizations are carried out until the total energy converges below 10-4 Ry and the atomic forces fall below 10-3Ry.
RESULTS AND DISCUSSION
Effect of K-Points Mesh Density on the Geometrical Structure of Graphene
The density of the k-points mesh in density functional theory (DFT) simulations plays a crucial role in determining the accuracy of structural parameters in periodic systems such as graphene 1. In this study, the k-points mesh was varied from 8×8×1 until 20×20×1, with and without the application of offset. The structural parameters analyzed include the lattice constants (lattice A and B) and the carbon–carbon (C–C) bond length.
For the 8×8×1 mesh, the lattice constants were found to be 2.4651 Å and 2.4643 Å with offset, and 2.4651 Å and 2.4648 Å without offset. The corresponding C–C bond lengths were 1.4231 Å and 1.4234 Å, respectively. At a higher mesh density of 12×12×1, a slight reduction in the C–C bond length was observed under offset conditions (1.4206 Å), accompanied by lattice constants of 2.4613 Å and 2.4665 Å. Under non-offset conditions, the structural parameters were more consistent, with lattice A = 2.4650 Å, lattice B = 2.4645 Å, and a C–C bond length of 1.4232 Å. At 16×16×1, the lattice constants became nearly identical (~2.4646 Å), with C–C bond lengths of 1.4240 Å (offset) and 1.4260 Å (no offset). At the highest mesh density (20×20×1), the structural parameters exhibited strong convergence, with lattice constants ~2.4648 Å and C–C bond lengths around 1.423 Å for both offset conditions.
In general, increasing the k-points mesh density leads to more stable and converged structural parameters. The differences between offset and non-offset conditions are relatively minor, with the maximum deviation in C–C bond length being less than 0.0054 Å. These slight discrepancies may arise from the sensitivity of Brillouin zone sampling to the initial k-point positions, particularly at lower mesh densities 14. However, at a mesh density of 16×16×1 and above, the structural parameters exhibit excellent consistency, indicating that structural convergence has been achieved. Therefore, the 16×16×1 mesh can be considered an optimal balance between computational cost and accuracy in this study.
Previous studies showed that, similar convergence test was considered reached around 6×6×1 (at 550 eV cutoff) and stabilized by 12×12×1 in practice 17,18. Meanwhile, in our findings (by using k-pints mesh 16×16×1 and above) align well, we see convergence of lattice (~2.4646Å) and bond length (~1.423Å), which matches or even improves upon those benchmarks. Consequently, we confirm that the 16×16×1 k‑point mesh delivers fully converged results that align with both theoretical and experimental benchmarks. Furthermore, these findings highlight the importance of selecting an adequately dense k-points mesh in DFT simulations to obtain reliable structural representations. This approach aligns with best practices in ab initio calculations, where numerical convergence is essential to ensuring the validity of simulation outcomes 19,20.
Total Energy Convergence of Graphene
Based on the obtained data, the total energy of the system was calculated using various k-points meshes, namely 8x8x1, 12x12x1, 16x16x1, and 20x20x1, under both offset and no-offset conditions. Generally, the total energy values ranged from -24.8623 to -24.8628 eV, exhibiting minimal variation across different k-points meshes and between offset conditions.
Figure 1 and 2 showed the visualization of total energy as a function of k-points mesh density indicates that significant energy convergence is achieved at the 12x12x1 mesh and remains stable up to 20x20x1, for both offset and no-offset cases. This finding suggests that increasing the k-points mesh density beyond 12x12x1 does not lead to substantial changes in total energy, thereby deeming this density sufficient for optimal computational accuracy. Energy convergence is critical as it reflects the stability of simulation results and prevents excessive computational resource usage without meaningful improvements in precision. A comparison between offset and no-offset conditions reveals nearly identical and stable total energy values, although the offset condition tends to require longer computational times. We note an increase in total energy at the 12×12×1 k-point mesh, which is consistent with the oscillatory convergence behavior documented in DFT literature 21. Discrete sampling of the Brillouin zone can alternately under- or over-estimate the integral over band energies, leading to non-monotonic energy trends. Such fluctuations are well known 21,22 and validated in high-throughput convergence frameworks 23.
Therefore, the choice between offset and no-offset conditions can be adjusted based on computational efficiency priorities without compromising energy accuracy. These results align with previous studies on energy convergence in DFT simulations using the SIESTA method 21, emphasizing the importance of selecting an appropriate k-points mesh to balance accuracy and computational efficiency.
CONCLUSION
This study has systematically evaluated the influence of k-points mesh density and offset conditions on the structural accuracy, total energy convergence, and computational efficiency of a pristine graphene system. The total energy results indicate that convergence is achieved at a k-points mesh of 12×12×1, with negligible variations up to 20×20×1. Similarly, structural parameters, including lattice constants and C–C bond lengths, demonstrate minimal deviation at higher mesh densities. However, computational time increases non-linearly with k-points density, especially under offset conditions, highlighting the trade-off between precision and computational cost. Based on our comprehensive assessment of energy stability, structural consistency, and time efficiency, the 16×16×1 no-offset k-points mesh emerges as the most balanced and reliable configuration. It yields the lowest total energy, exhibits excellent agreement with established structural benchmarks, and avoids excessive computational demand. This makes it particularly suitable as a reference system for future ab initio studies, such as H2O adsorption on graphene, where accurate baseline energies are critical for computing adsorption energetics. The findings underscore the importance of k-points convergence testing in density functional theory (DFT) simulations and support prior literature emphasizing the balance between computational accuracy and efficiency. Future adsorption studies can confidently adopt the identified k-points mesh to ensure both reliable results and computational feasibility.