Investigate switching integrator per particle

[RudolfWeeber]

Right now, the integrator is chosen globally (via integ_switch). Specific intregrators contain their own loop over all particles.

This has two weaknesses:

• some common functions like Verlet list updates have to copied into each integrator, or an extra iteration over all particles is incurred, if Verlet list management is to be separated
• Inertialess tracers are currently implemented as virtual sites, although, they conceptually are more like an integrator (#2695, #3059).

Using definitions:

• virtual sites ar particles who's properties position/orientation/velocity are dreived from other particles
• integrators implement rules that propagate the particles

We should investigate, whether

• exposing integrator per-particle kernels rather than sweeps
• putting the sweep over all particle in the central integration function and using switching kernels on a per-particle level
• Verlet list criterion check can then be integrated as an additional per-particle kernel
• particles can have an additional flag which overrides the default integrator This would be used for the inertialess tracers

In a second step (or in one go), integrator parameters can be put into classes. This would be time step for velocity verlet, pre-factor and cap for steepest descent, etc.

@fweik, this is how I understood our discussions on the topic. Is there something missing?

[fweik]

Yes that sounds reasonable. If you want to coerce the virtual sites into that abstraction (they are really not per particle), than they will always need special treatment. Also it has to be considered if rattle can be integrated into this framework.

For the virtual sites I think abstractly what you want to do is in the first sweep defer some particles for later integration. You could e.g. have integrators return either new positions for a particle or a defer tag, in which the particle is noted for later propagation. Then, if there are deferred particles, you do the next loop only with the left ones and so on until there are none left. Potentially rattle also fits into this construction e.g. rattle could defer particles until it is happy with the tolerance.

[jngrad]

Equations of motion for integrators and thermostats

Used notation:

• position: $x(t)$
• quaternion: $q(t)$
• velocity: $v(t)$
• angular velocity: $\omega(t)$
• force: $f(t)$
• torque: $\tau(t)$
• mass: $m$
• moment of inertia: $I$
• interpolated fluid velocity: $u(x(t))$
• box length: $L(t)$
• time step: $dt$
• temperature: $T$
• friction: $\gamma$
• noise: $\eta()$
• quaternion rotation matrix: $R()$

Velocity Verlet

• documentation
• step 1: requires time_step, skin (to resort particles) $\omega(t + dt / 2) = \omega(t) + \tau(t) / I \cdot dt / 2 + \Delta\omega \cdot dt / 2$ $v(t + dt / 2) = v(t) + f(t) / m \cdot dt / 2$ $x(t + dt / 2) = x(t) + v(t + dt / 2) \cdot dt$ $q(t + dt / 2) = q(t) + \Delta q \cdot dt + \Delta \Delta q \cdot dt^2 / 2 - \lambda q(t)$
• step 2: requires time_step $\omega(t + dt) = \omega(t + dt / 2) + \tau(t + dt) / I \cdot dt / 2 + \Delta\omega \cdot dt / 2$ $v(t + dt) = v(t + dt / 2) + f(t + dt) / m \cdot dt / 2$ $x(t + dt) = x(t + dt / 2)$ $q(t + dt) = q(t + dt / 2)$

Velocity Verlet NpT

• documentation
• frictionless rotation
• step 1: requires time_step, temperature, npt_iso, box_length $\omega(t + dt / 2) = \omega(t) + \tau(t) / I \cdot dt / 2 + \Delta\omega \cdot dt / 2$ $v(t + dt / 2) = v(t) + f(t) / m \cdot dt / 2 + (\gamma \cdot v(t) \cdot dt / 2 + \sigma(T, \gamma, dt) \cdot \eta()) / m$ $x(t + dt / 2) = \frac{L(t + dt)}{L(t)} \cdot (x(t) + \frac{L^2(t)}{L^2(t + dt / 2)} \cdot v(t + dt / 2) \cdot dt )$ $v(t + dt / 2) = \frac{L(t)}{L(t + dt)} \cdot v(t + dt / 2)$ $q(t + dt / 2) = q(t) + \Delta q \cdot dt + \Delta \Delta q \cdot dt^2 / 2 - \lambda q(t)$
• step 2: requires time_step, npt_iso $\omega(t + dt) = \omega(t + dt / 2) + \tau(t + dt) / I \cdot dt / 2 + \Delta\omega \cdot dt / 2$ $v(t + dt) = v(t + dt / 2) + f(t + dt) / m \cdot dt / 2 + (\gamma \cdot v(t + dt / 2) \cdot dt / 2 + \sigma(T, \gamma, dt) \cdot \eta()) / m$ $x(t + dt) = x(t + dt / 2)$ $q(t + dt) = q(t + dt / 2)$

Stokesian Dynamics

• documentation
• step 1: requires time_step, skin (to resort particles), sd_viscosity, sd_kT, sd_flags, radius_dict $x(t + dt) = x(t) + v(t + dt) \cdot dt$ $q(t + dt) = q(t) + R(\omega(t + dt) \cdot dt)$
• step 2: nothing

Brownian Dynamics

• documentation
• step 1: requires time_step, temperature, skin (to resort particles) $x(t + dt) = x(t) + f(t) / \gamma \cdot dt$ $q(t + dt) = q(t) + R(\tau(t) / \gamma \cdot dt)$ $v(t + dt) = f(t) / \gamma$ $\omega(t + dt) = \tau(t) / \gamma$ $x(t + dt) = x(t + dt) + \sigma(T, \gamma, dt) \cdot \eta() \cdot dt^{1/2}$ $q(t + dt) = q(t + dt) + R(\sigma(T, \gamma, dt) \cdot \eta() \cdot dt^{1/2})$ $v(t + dt) = v(t + dt) + \sigma(T, m) \cdot \eta() \cdot m^{1/2}$ $\omega(t + dt) = \omega(t + dt) + \sigma(T, m) \cdot \eta() \cdot I^{1/2}$
• step 2: nothing

Steepest descent

• documentation
• step 1: requires params $x(t + dt) = x(t) + \min(\gamma f(t), r^{\mathrm{max}})$ $q(t + dt) = q(t) + R(\min(\gamma \tau(t), r^{\mathrm{max}}))$
• step 2: nothing

Langevin

• documentation
• forces: requires time_step, temperature, langevin $f(t) = f^{\mathrm{ext}}(t) - \gamma v(t) + \eta(T, \gamma, dt)$ $\tau(t) = \tau^{\mathrm{ext}}(t) - \gamma \omega(t) + \eta(T, \gamma, dt)$

LB

• documentation
• forces: requires time_step, temperature, gamma_lb $f(t) = f^{\mathrm{ext}}(t) - \gamma \cdot (v(t) - u(x(t))) + \sigma(T, \gamma, dt) \cdot \eta()$

[RudolfWeeber]

@jngrad, thanks for collecting the status quo.

In the future, there should be

• possibility to setup several propagation schemes for translation and rotation. For the most part, schemees should work either on translation or rotation.
• Possibility to choose one default propagation for translation and one for rotation. This applies to all particles, which do not have an other selection.
• Possibility to select for specific particles one of the defined propagation schemes for translation and one for rotation

There will likely be the following propagation schemes:

• Translation
  • Velocity Verlet (dt, particle specific mass)
  • Velocity Verlet + Langevin (dt, kT, gamma, particle specific mass)
  • Velocity Verlet + LB Langevin coupling (dt, kT, gamma, particle specific gamma, mass)
  • Brownian dynamics, (dt, gamma, particle specific gamma)
  • Stokesian dynamics (dt, gamma, eta, radii, particle speciffic gamma) Don't know if this also acts on rotation @biermanncarl?
  • Velocity verlet + NPT box size rescaling (dt, gamma, gamma_box, piston_mass, mass on the particle)
  • Virtual sites relative (none globally, base particle id distance and quaternion on the particles)
  • Inertialess LB tracers (IBM method) (none probably) @sgekle
  • Virtual sites at center of triangle (replaces removed oif out direction) (None globally, base particle ids on the particles) @icimrak
  • Not sure how rattle fits in. Probably only as an addon to Velocity Verlet
  • ignore
• Rotation
  • Velocity Verlet (dt, particle specific inertia)
  • Velocity Verlet + Langevin (dt, kT, gamma, particle specific gamma, inertia)
  • Brownian dynamics (dt, gamma, particle specific gamma)
  • VS relative orientation update (no globals, id of base particle, rel orientation quaternion)
  • ignore

[RudolfWeeber]

• In some cases VS are thermalized (e.g. raspberries)
• Some propagators are not single particle functions (Stokesian)

[RudolfWeeber]

Single particle propagators without pre/post tasks

Below propagators have kernels of the form

new_pos, new_vel = kernel(particle)

force is optional in some cases. Analogous for rotation. There is (to my understanding) no task which has to run before they can be applied

• Translation
  • Velocity Verlet (dt, particle specific mass)
  • Velocity Verlet + Langevin (dt, kT, gamma, particle specific mass)
  • Velocity Verlet + LB Langevin coupling (dt, kT, gamma, particle specific gamma, mass). Adds LB forces as side effect
  • Brownian dynamics, (dt, gamma, particle specific gamma)
• Rotation
  • Velocity Verlet (dt, particle specific inertia)
  • Velocity Verlet + Langevin (dt, kT, gamma, particle specific gamma, inertia)
  • Brownian dynamics (dt, gamma, particle specific gamma)

Single particle propagators with pre/post tasks

The following propagators have single particle kernels + a pre- and/or a post-task

• Velocity verlet + NPT box size rescaling (dt, gamma, gamma_box, piston_mass, mass on the particle)
  • The propagator also produces a new box size
  • All virtual sites (other particles)
  • Requires the base particles to be already propagated.
  • Inertialess LB tracers (IBM method) (adds lb forces as side effect)
  • Require ghost force reduction

Collective propagators

• Stokesian dynamics (dt, gamma, eta, radii, particle speciffic gamma)
  • Don't know if this also acts on rotation @biermanncarl?

Not sure how rattle fits in. Probably only as an addon to Velocity Verlet

[RudolfWeeber]

So, if we only had single particle kernels which work independently, we'd have

for prop in active_propagatros:
    prop.pre(filter(particles,prop))
for part in particles:
    select_propagator(part).propagate(part) 
for prop in propagatros:
    prop.post(filter(particles,prop))

Here, filter() would filter the particles that qualify for a specific propagator. select_propagatro() will either reutrn the default (system-wide) propagator or the one specified on the single particle.

Alternatively, one can have one traversal of particles per propagator

for prop in all_propagators:
    prop.propagate(filter(particles,prop))

The latter approach works also with collective particles and would allow us to specify an order (e.g., virtual sites last) The single particle propagators could then share (part of) the propagate() implementation.

The first approach would allow to combine the vve rlet list update check into the propagation traversal. Otherwise, an independent pass is needed. (Currently, verlet list checks are copied into all propagatros, which I don't think, we want to keep).

We also need to decide, how functional we want this to be. I.e., which of the following should the propagator change "in-place"and which whoudl be returned as some state-change datastructure, which is acted upon later?

• position, quaternion
• velocity, torque
• box size (hello npt)
• forces on lb