Navigation :

Detailled example: Verlet

Work in progress!

The Verlet algorithm

According to Wikipedia, the Verlet algorithm is defined as:

Calculate $\vec{x}(t + \Delta t) = \vec{x}(t) + \vec{v}(t) \Delta t + \tfrac12 \vec{a}(t) \Delta t^2$ .
Derive $\vec{a}(t + \Delta t)$ from the interaction potential using $\vec{x}(t + \Delta t)$ .
Calculate $\vec{v}(t + \Delta t) = \vec{v}(t) + \tfrac12 \big(\vec{a}(t) + \vec{a}(t + \Delta t)\big)\Delta t$ .

The Verlet operator is represented by the type propagator_verlet_t which extends propagator_t:

Definition of propagator_verlet_t

Note, that the class definition does not add anything to the propagator_t class. The only differences are the definition of the operations, and the overloaded constructor.

For the Verlet propagator, we need to define the following operations:


  ! Specific verlet propagation operations identifiers
  character(len=30), public, parameter ::      &
    VERLET_START       = 'VERLET_START',       &
    VERLET_FINISH      = 'VERLET_FINISH',      &
    VERLET_UPDATE_POS  = 'VERLET_UPDATE_POS',  &
    VERLET_COMPUTE_ACC = 'VERLET_COMPUTE_ACC', &
    VERLET_COMPUTE_VEL = 'VERLET_COMPUTE_VEL'

  ! Specific verlet propagation operations
  type(algorithmic_operation_t), public, parameter :: &
    OP_VERLET_START       = algorithmic_operation_t(VERLET_START,       'Starting Verlet propagation'),               &
    OP_VERLET_FINISH      = algorithmic_operation_t(VERLET_FINISH,      'Finishing Verlet propagation'),              &
    OP_VERLET_UPDATE_POS  = algorithmic_operation_t(VERLET_UPDATE_POS,  'Propagation step - Updating positions'),     &
    OP_VERLET_COMPUTE_ACC = algorithmic_operation_t(VERLET_COMPUTE_ACC, 'Propagation step - Computing acceleration'), &
    OP_VERLET_COMPUTE_VEL = algorithmic_operation_t(VERLET_COMPUTE_VEL, 'Propagation step - Computing velocity')

These are used to define the algorithm, which is done in the constructor of the propagator:

  function propagator_verlet_constructor(dt) result(this)
    FLOAT,                     intent(in) :: dt
    type(propagator_verlet_t), pointer    :: this

    PUSH_SUB(propagator_verlet_constructor)

    allocate(this)

    this%start_operation = OP_VERLET_START
    this%final_operation = OP_VERLET_FINISH

    call this%add_operation(OP_VERLET_UPDATE_POS)
    call this%add_operation(OP_UPDATE_COUPLINGS)
    call this%add_operation(OP_UPDATE_INTERACTIONS)
    call this%add_operation(OP_VERLET_COMPUTE_ACC)
    call this%add_operation(OP_VERLET_COMPUTE_VEL)
    call this%add_operation(OP_ITERATION_DONE)
    call this%add_operation(OP_REWIND_ALGORITHM)

    ! Verlet has only one algorithmic step
    this%algo_steps = 1

    this%dt = dt

    POP_SUB(propagator_verlet_constructor)
  end function propagator_verlet_constructor

The algorithm also uses steps, which are not specific to the Verlet algorithm, and are defined in the propagator_oct_m module.

Note, the difference between OP_VERLET_FINISH and OP_FINISHED. The former denotes a specific step, to be taken at the end of one time step, while the latter generally denotes the end of a time step.

As can be seen in this example, the definition of the propagator in terms of the algorithmic operations is a direct translation of the algorithm, shown above.

Implementation of the steps

So far, the propagator is only defined in an abstract way. It is down to the individual systems (or better, their programmers) to implement the individual steps of the algorithm, in terms of the dynamic variables for that system. An easy examample to demonstrate this is the classical particle, as implemented in classical_particles_t

definition of classical_particles_t

This describes an extremely simple system, consisting of a set of classical, neutral particles. The dynamic variables (i.e. the state) of the system are the positions, velocities and accelerations.

One ‘‘tick’’ of the propagator is defined in the function classical_particle_do_td. As can be seen in the code below, this function implements all possible algorithmic steps for all propagators, allowed for that system.

Example implementation for classical particles

The timeline explained

The Verlet algorithm is a good (because simple) example to illustrate how the systems and the interaction are updated as the code progresses through the main loop.

0: multisystem_propagation_start

System: root

Clocks:	In	Out
system:	0.00000E+00	0.00000E+00
prop:	0.00000E+00	0.00000E+00

This graph illustrates how the state machine is stepping through the algorithm. Each system is picking the next algorithmic step from the propagator. For the containers (i.e. ‘‘root’’ and ‘‘earth’'), the only steps are ‘‘Updating interactions’’ and ‘‘Finished’’. The real systems, on the other hand, are progressing orderly through the operations, defined in the propagator.

Example with different time steps