# softmatter:Velocity Verlet Algorithm

The original Verlet Algorithm does not use the velocity to compute the new position; therefore the velocity is only derived from the old and new positions. Velocity Verlet algorithm utilizes the new forces to update velocities at time (t + Δt), allowing the position and velocity are automatically synchronized. This velocity Verlet algorithm proceeds as follows 1:

## General Algorithm

Expand the position of a particle at time t in a truncated Taylor series:

$\mathbf{r}(t + \Delta t) = \mathbf{r}(t) + \mathbf{v}(t) \Delta t + \frac{\mathbf{f}(t)}{2m} (\Delta t)^2 \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad [1]$

where $m, \mathbf{r}(t), \mathbf{v}(t)$ are the mass, position and velocity of the particle, $\mathbf{f}(t)$ is the force acting on the particle at time t and Δt is the time step of integration.

Unlike the original Verlet Algorithm, the velocity now is updated from the forces at the previous and current time steps, i.e.

$\mathbf{v}(t + \Delta t) = \mathbf{v}(t) + \frac{\mathbf{f}(t + \Delta t) + \mathbf{f}(t)}{2m} (\Delta t) \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad [2]$

This integration scheme is indeed equivalent to the original scheme [1] when we write down the expression of position at the next time step (t + 2Δt)

$\mathbf{r}(t + 2 \Delta t) = \mathbf{r}(t + \Delta t) + \mathbf{v}(t + \Delta t) \Delta t + \frac{\mathbf{f}(t + \Delta t)}{2m} (\Delta t)^2 \qquad \qquad \qquad \qquad \qquad \qquad [3]$

and re-arrange equation [1] into

$\mathbf{r}(t) = \mathbf{r}(t + \Delta t) - \mathbf{v}(t) \Delta t - \frac{\mathbf{f}(t)}{2m} (\Delta t)^2 \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad [4]$

By addition of [3] and [4], we get

$\mathbf{r}(t + 2 \Delta t) + \mathbf{r}(t) = 2\mathbf{r}(t + \Delta t) + [\mathbf{v}(t + \Delta t) - \mathbf{v}(t)]\Delta t + \frac{\mathbf{f}(t + \Delta t) - \mathbf{f}(t)}{2m} (\Delta t)^2 \qquad \qquad [5]$

Finally, substitution of equation [2] for $\mathbf{v}(t + \Delta t)$ into [5] gives us the coordinate equation of the Verlet Algorithm:

$\mathbf{r}(t + 2 \Delta t) + \mathbf{r}(t) = 2\mathbf{r}(t + \Delta t) + \frac{\mathbf{f}(t + \Delta t)}{2m} (\Delta t)^2$

This result shows that velocity Verlet algorithm yields the same order of accuracy in velocity as in the original scheme, i.e. t)2. Higher order of velocity accuracy can be achieved when additional corrections are added, e.g Beeman algorithm or Velocity-corrected Verlet Algorithm.

### Example Source Code

The velocity Verlet algorithm is implemented based on equation [1] and [2]. Given initial positions and velocities, the trajectories of particles will be determined as the simulation progresses. Additionally initial force should be also computed to boot up the integrating loop. A typical integration routine is given below 2:

void Integrate()
{
// Initialize positions and velocities
Initialize();

// Compute the forces from the initial positions and velocities
ComputeForce();

for (int timestep = _begin; timestep <= _end; timestep++)
{
// Initial integration using the forces from the previous step, noting that velocities are advanced in half time step
for (int i = 0; i < Nparticles; i++)
{
v[i] = v[i] + (f[i] / mass[i]) * delta_t / 2;
x[i] = x[i] + v[i] * delta_t;
}

// Compute new forces from new positions
ComputeForce();

// Final integration, update the velocities by half time step using the new forces
for (int i = 0; i < Nparticles; i++)
v[i] = v[i] + (f[i] / mass[i]) * delta_t / 2;
}

}


• In this implementation, the update of velocities is divided into 2 steps. The initial integration uses the forces from the previous time step to advance velocity by half time step

$\mathbf{v}(t + \Delta t / 2) = \mathbf{v}(t) + \frac{\mathbf{f}(t)}{2m} (\Delta t)$

The new position is computed using equation [1] where the last two terms in the right hand side is exactly the half time step velocity $\mathbf{v}(t + \Delta t / 2)$ obtained above

$\mathbf{r}(t + \Delta t) = \mathbf{r}(t) + \mathbf{v}(t) \Delta t + \frac{\mathbf{f}(t)}{2m} (\Delta t)^2 = \mathbf{r}(t) + \mathbf{v}(t + \Delta t / 2) (\Delta t)$

The final integration completes the velocity update in equation [2], after the new forces are computed:

$\mathbf{v}(t + \Delta t) = \mathbf{v}(t + \Delta t / 2) + \frac{\mathbf{f}(t + \Delta t)}{2m} (\Delta t)$

## Algorithm for Rigid Bodies

The velocity Verlet algorithm can be extended to the motion of rigid particles (for more on rigid body motion, see quaternions). For rigid object motion, the rotational degrees of freedom are governed by

$\frac{d\mathbf{L}(t)}{dt} = \frac{d(\mathbf{I}(t) \omega(t))}{dt} = \mathbf{T}(t)$

where $\mathbf{L}(t)$, $\mathbf{\omega}(t)$ and $\mathbf{I}(t)$ are the angular momentum, angular velocity and moment of inertia tentor of the rigid body, respectively. $\mathbf{T}(t)$ is the torque applied to the rigid body. All of these quantities are in space-fixed coordinates. The angular momentum can be updated in the same manner as translational velocity. For more details, please see Rigid Body Motion.

The following listing illustrates a typical implementation of velocity Verlet algorithm for rigid body motion 2.

### Example Source Code

void Integrate_Rigidbody()
{
// Initialize positions and velocities
Initialize();

// Compute the forces from the initial positions and velocities
ComputeForceAndTorque();

for (int timestep = _begin; timestep <= _end; timestep++)
{
// Initial integration using the forces from the previous step, noting that velocities are advanced in half time step
for (int i = 0; i < Nbodies; i++)
{
// Translational part
vm[i] = vm[i] + (fm[i] / mass[i]) * delta_t / 2;
xm[i] = xm[i] + vm[i] * delta_t;

// Rotational part
angular_momentum[i] = angular_momentum[i] + torque[i] * delta_t / 2;

// Update the position, linear and angular velocities of body i
SetPositionAndVelocity(i);
}

// Compute new forces from new positions
ComputeForceAndTorque();

// Final integration, update the velocities by half time step using the new forces
for (int i = 0; i < Nbodies; i++)
{
// Translational part
vm[i] = vm[i] + (fm[i] / mass[i]) * delta_t / 2;

// Rotational part
angular_momentum[i] = angular_momentum[i] + torque[i] * delta_t / 2;

// Update the linear and angular velocities of body i
SetVelocity(i);
}
}
}