Wavefunction extrapolation greatly reduces the number of self-consistent field(SCF)iterations and thus the overall computational cost of Born-Oppenheimer molecular dynamics(BOMD)that is based on the Kohn–Sham density functional theory.Going against the intuition that the higher order of extrapolation possesses a better accuracy,we demonstrate,from both theoretical and numerical perspectives,that the extrapolation accuracy firstly increases and then decreases with respect to the order,and an optimal extrapolation order in terms of minimal number of SCF iterations always exists.We also prove that the optimal order tends to be larger when using larger MD time steps or more strict SCF convergence criteria.By example BOMD simulations of a solid copper system,we show that the optimal extrapolation order covers a broad range when varying the MD time step or the SCF convergence criterion.Therefore,we suggest the necessity for BOMD simulation packages to open the user interface and to provide more choices on the extrapolation order.Another factor that may influence the extrapolation accuracy is the alignment scheme that eliminates the discontinuity in the wavefunctions with respect to the atomic or cell variables.We prove the equivalence between the two existing schemes,thus the implementation of either of them does not lead to essential difference in the extrapolation accuracy.