Below, we decided to calculate the velocities (when required) directly. This changed the to do list for this issue: - [x] Derive velocity expressions and update the documentation - [x] Implement new velocity calculation, to be performed only if trajectory writing is about to happen - [x] Remove previous workarounds needed to keep positions stable during update - [x] Move virtual site calculation to top of the loop right before force calculation, avoiding a) problems with pressure scaling, b) reducing the likelihood of interference with other algorithms between the moment they are calculated and the moment they enter the dynamic properties (forces), and c) remove the need for any initialization before the loop. - [x] Introduce tests (see below) - [x] Check compatibility with GPU integration - [ ] Check possible bug in MiMiC I'm leaving the original description below for reference --- Virtual sites are massless "particles" whose positions are fully determined by a function of other (non-virtual or virtual) atoms. They may be involved in interactions, but forces acting on them are redistributed to particles having mass. Their positions are updated after non-virtual particles were updated and constrained. Their velocities are not used in their propagation, but are calculated from their displacement. This suggests that * virtual site positions should only be changed by the virtual site algorithm, * virtual site velocities are zero at the first step, and calculated from the displacement of the positions after that, and should only be changed by the virtual site algorithm. Calculation of the velocity from the position displacement is done as `v = (newX - oldX) / dt`. This results in half-step velocities as expected in leap-frog. Velocity verlet integrator reports atom velocities at full steps, which can't be calculated for the virtual sites until after the next position update. This would be an argument not to report the virtual site velocities. Current md-vv is setting the velocities to zero shortly before trajectory writing, but it's unclear whether this is by design or by accident. Reporting zero velocities can be confusing to users, however - the virtual site positions might seem frozen. - [ ] **Action item 1.1: Decide if virtual site velocities should be reported at all - almost all issues below could easily be solved by simply zero-ing the virtual velocities right before trajectory output... What is the use case for virtual site velocities?** - [ ] **Action item 1.2: Decide if md-vv virtual site velocities should be half-step velocities or explicitly zero-ed.** - [ ] **Action item 1.2a: Adapt implementation to decision of 1.2.** - [ ] **Action item 1.3: Implement test checking that virtual site velocities conform to expectation.** Calculating the velocity from the displacement requires that the **virtual positions are not updated outside of the virtual site algorithm** (e.g. during the propagation of regular atoms). To avoid that virtual site positions get updated, algorithms can either check each atom for its particle type and skip propagation (done in md-vv, sd, bd, which also zero the velocities - probably not necessary) or zero the virtual site velocities before update (done by the simple and simple SIMD versions of leap-frog, such that the update [now with zero velocities and zero forces] can be vectorized). The general leap-frog propagation does neither of the two, resulting in wrong virtual site velocities. - [ ] **Action item 2: Fix the general leap-frog update to skip virtual sites.** The virtual site positions and velocities are recalculated after the regular propagation, but they are changed again in `update_pcouple_after_coordinates` when using Berendsen (positions only), C-Rescale (positions and velocities) and Parrinello-Rahman (positions only). Scaling the virtual site positions as the normal particle positions can lead to **small errors in the virtual site positions**, as not all virtual sites are calculated as linear combinations of other particles. - [ ] **Action item 3: Move virtual site calculation after pressure coupling (make sure positions are final before updating virtual sites).** Even after moving the virtual site calculation after the pressure coupling scales the coordinates, having the pressure coupling scaling the virtual site positions will cause the virtual site velocities not to be exactly `v = (newX - oldX) / dt`. - [ ] **Action item 4: Exclude virtual site positions from pressure coupling position scaling.** The virtual site velocities are scaled during COM motion removal. Since the velocities represent the previous propagation, and virtual sites don't have a mass, they should probably not be affected by COM motion removal. - [ ] **Action item 5: Exclude virtual site velocities from COM motion removal.** The virtual site velocities are also rescaled during temperature coupling for md-vv. For non-Trotter schemes, this does not matter because it happens right before coordinate propagation (which zeroes the velocities), but has some potential for future bugs. For Trotter, this might influence the velocities if we decide to report them (see action item 1). - [ ] **Action item 6.1: Exclude virtual site velocities from md-vv Trotter temperature coupling rescaling, unless they are defined to be zero (Action item 1.2).** - [ ] **Action item 6.2: Decide if we need to exclude virtual site velocities from temperature coupling rescaling in other cases.** Subtask of #3849.