softmatter:Velocity Verlet Algorithm

From NSDL Materials Digital Library Soft Matter Wiki

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
- 1.1 Example Source Code
2 Algorithm for Rigid Bodies
- 2.1 Example Source Code
3 References

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 ²:

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;
     }

}

Comments:

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 ².

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); 
          }
     }
}

Comments

SetPositionAndVelocity(i) updates the position and velocity of constituent components of the rigid body i and its angular velocity.

SetVelocity(i) updates the velocity of constituent components of the rigid body i and its angular velocity.

References

[1] Understanding molecular simulation : from algorithms to applications
[2] LAMMPS source code, LAMMPS website

Retrieved from "http://matdl.org/matdlwiki/index.php/softmatter:Velocity_Verlet_Algorithm"

Categories: Softmatter Simulation Methods | Softmatter Integration Schemes

softmatter:Velocity Verlet Algorithm

From NSDL Materials Digital Library Soft Matter Wiki

Contents

General Algorithm

Example Source Code

Algorithm for Rigid Bodies

Example Source Code

References

Views

Personal tools

MatDL Wiki

Search

Toolbox