The finite difference method can generate MD trajectories with continuous potential models. The essential idea is that the integration is divided into many small steps, each separated by a fixed time interval . The total interactions on each particle at time can be calculated from the sum of interactions from other particles. The force is assumed to be constant during the time step and . The forces and accelerations of the particles in new positions can be determined, and so on.

There are many algorithms for integrating the equations of motion using finite difference methods, which are commonly used in MD simulations. All algorithms assume that the positions and dynamics properties (velocities, accelerations, etc.) can be approximated in Taylor series expansions:

The execution of the Verlet algorithm is straightforward and comprises , , and the accelerations . One of its disadvantage is that the positions are obtained by adding a small term, , to the difference of two much larger terms, and ). This may lead to a loss of precision. The lack of an explicit velocity term in the equations makes it difficult to obtain the velocities. The velocities are not available until the positions have been computed at the next step. In addition, it is not a self-starting algorithm; the new positions are obtained from the current positions and the previous position . At there is obviously only one set of positions and so it is necessary to employ some other means to obtain positions at . One way to obtain is to use the Taylor series, Eq. (9), truncated after the first term. Thus, .

The velocity Verlet method gives positions, velocities and accelerations at the same time and does not influence the precision:

The velocity Verlet method [26] is actually implemented as a three-stage procedure because, as can be seen from Eq. (18), calculating the new velocities requires the accelerations at both and . Thus in the first step the positions at are calculated according to Eq. (17) using the velocities and the accelerations at time . The velocities at time are then determined using