The phenomena of Coulomb explosion require the consideration of special relativity due to the involvement of high energy electrons or ions. It is known that laser ablation processes at high laser intensities may lead to the Coulomb explosion, and their released energy is in the regime of kEV to MeV. In contrast to conventional molecular dynamics (MD) simulations, we adopt the three-dimensional relativistic molecular dynamics (RMD) method to consider the effects of special relativity in the conventional MD simulation for charged particles in strong electromagnetic fields. Furthermore, we develop a Coulomb force scheme, combined with the Lennard-Jones potential, to calculate interactions between charged particles, and adopt a Verlet list scheme to compute the interactions between each particle. The energy transfer from the laser pulses to the solid surface is not directly simulated. Instead, we directly assign ion charges to the surface atoms that are illuminated by the laser. By introducing the Coulomb potential into the Lennard-Jones potential, we are able to mimic the laser energy being dumped into the xenon (Xe) solid, and track the motion of each Xe atom. In other words, the laser intensity is simulated by using the repulsive forces from the Coulomb potential. Both nonrelativistic and relativistic simulations are performed, and the RMD method provides more realistic results, in particular, when high-intensity laser is used. In addition, it is found that the damage depth does not increase with repeated laser ablation when the pulse frequency is comparable to the duration of the pulse. Furthermore, we report the time evolution of energy propagation in space in the laser ablation process. The temporal-spatial distribution of energy indirectly indicates the temperature evolution on the surface of the Xe solid under intense laser illumination. [ABSTRACT FROM AUTHOR]