Dear X.J.Xu
Such total energy drifts indicate that the results of the simulation are not reliable, the errors in the integration of the equation of motion are to big.
The conserved quantity must be conserved.
The increase in time per MD step is most probably related to an increase in the number of SCF iterations, if you are running BOMD.
If this is the case, i.e., if it becomes more and more difficult to optimise the wave function, it might indicate that the computed forces are rather wrong and lead to the sampling of weird geometries.
By the way, have you considered that graphene is a semi-metal? This makes the convergence of the wave function even more problematic.
Before worrying on how to speed up the sampling, it is recommendable to find settings that produce reliable results, for instance further reducing EPS_SCF, adding virtual MOS, increasing the cutoff, etc.
Kind regards
Marcella