How energetic particles (EPs) affect the internal kink mode (IKM) is a critical question for the stable operation of future burning plasma devices like ITER and the EU DEMO. In this paper, the IKM in tokamak geometry is numerically systematically investigated utilizing the full toroidal magnetohydrodynamic (MHD) code MARS-F and the non-perturbative MHD-kinetic hybrid code MARS-K. Emphasis is placed on the role of kinetic effects arising from the precessional drift motion of trapped EPs. Multiple plasma parameters are scanned including the ideal wall position d / a, the on-axis pressure ratio \hat\beta_0^* and the energy ratio \alpha_\mathrmT between the energetic and thermal particles, as well as the radial profiles of EP density and pressure. The kinetic effect is found to be either stabilizing or destabilizing depending on the EP equilibrium profiles. With the latter being constant fractions of the corresponding thermal particle profiles, the EP kinetic effect mildly destabilizes the IKM at low d / a and \alpha_\mathrmT. The mode growth rate monotonically increases with \beta_0^*, peaking at maximum d / a and \alpha_\mathrmT for \beta_0^*=0.5. For radially distributed profiles for the EP density and pressure fractions, however, the effect is stabilizing throughout the parameter spaces in \left(d / a, \alpha_\mathrmT\right). The enhanced stabilizing effect comes from the more peaked radial profiles for the EP density and pressure, which increases the magnitude of the resonance operator-given in Eq. (8)-by increasing the diamagnetic frequency of the EPs. At large \alpha_\mathrmT, the mode growth rate non-monotonically depends on \beta_0^*, reaching the peak value at maximum d / a but minimum \alpha_\mathrmT for \beta_0^*=0.5. This non-monotonic behavior stems from the fact that a higher mode frequency is excited by precessional drift resonance of trapped EPs for radially distributed profiles, while the mode frequency appears in both the numerator and denominator of the resonance operator. Significant quantitative differences are also computed, over the entire parameter space, between the aforementioned two profile models for EPs. Comparative analysis shows approximately one order of magnitude difference in the computed mode frequency, and even sign reversal, between the two EP profile models with increasing \alpha_\mathrmT. With the radial profiles of EP density and pressure more representative of the experimentally conditions, the mode frequency excited by EPs reaches an amplitude on the order of 10^-2, corresponding to real frequencies of several to tens of kHz . This result aligns well with experimental and numerical observations from multiple tokamaks, including DIII-D, EAST, and HL-2A. The amplitude increases monotonically with EP energy, further supporting the consistency with these observations. In experiments, the radial profiles of EP pressure and density can serve as key control parameters. By tailoring the deposition profile of auxiliary heating, the radial distribution of EPs can be actively manipulated, thereby enabling effective control of sawtooth oscillations. These results enhance understanding of the EP kinetic effects on IKM and provide physics insight for analyzing associated instabilities in future fusion devices.